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SUMMARY 

This report describes work performed on the development of a hierarchical 
real-time algorithm for optimal three-dimensional aircraft maneuvers using 
Singular Perturbation Theory (SPT). New theoretical results justify and develop 
systematic methods for real-time computation of nonlinear feedback controls by 
means of SPT and provide an assessment of the accuracy of the resulting SPT 
control. Practical results apply SPT to obtain a real-time feedback law for 
the three-dimensional minimum time long range intercept problem for an F-4 air- 
craft model (six state, three control variable, point mass model). Nonlinear 
feedback laws are presented for computing the optimal control variables u (throttle), 
a (bank angle) and a (angle-of-attack) as a function of target and pursuer air- 
craft states and desired terminal conditions. A real-time capabil ity assess- 
ment of the SPT algorithm on a TI9900 microcomputer has been performed and the 
control update rates have been determined. The storage and computational re- 
quirements of the algorithm are found to be well suited for on-board real-time 
implementation on a microcomputer. 

The accuracy of the SPT solution is analyzed and it is shown how "continua- 
tion-type" methods may be used to obtain exact optimal trajectories starting 
from the SPT solution. The advantage of using predictive terms to supplement 
the SPT feedback laws is demonstrated for the aircraft trajectory optimization 
problem. In particular, it is shown that the SPT approximation breaks down near 
the terminal phase and must be corrected by "continuation" and Generalized Mul- 
tiple Scale (GMS) methods. 


CHAPTER 1 
INTRODUCTION 


1*1 Summary of Contents 

This report describes work performed on the development of a real-time 
algorithm for optimal three-dimensional aircraft maneuvers using Singular 
Perturbation Theory (SPT). The optimization problem considered is that of 
minimum time long range interception. Nonlinear feedback laws are presented 
for computing the optimal control variables u (throttle), a (bank angle) and 
a (angle-of-attack) as a function of target and pursuer aircraft states. A 
real-time capability assessment of the SPT algorithm on a TI9900 microcomputer 
has been performed and the control update rates have been determined.' The 
storage and computational requirements of the algorithm are found to be well 
suited for on-board real-time implementation on a TI9900 microcomputer. 

_The organization of the report is as follows. Chapter 2 presents the air- 
craft model (a six-state, point mass approximation) considered in the project. 

The dynamic equations, constraints and numerical simulation are discussed in 
detail. Chapter 3 formulates the optimization problem and discusses exact methods 
of solution — in particular, the continuation method. Chapter 4 presents the 
general theory of feedback control law computation using SPT. The problems of 
computational efficiency and accuracy of the SPT method are discussed in detail. 
Chapter 5 applies these theoretical results to the aircraft optimization problem 
formulated in Chapter 3. We detail the computational procedures and the neces- 
sary approximations made to obtain computationally feasible solutions. Chapter 
6 describes the final real-time algorithm and presents numerical examples for 
various flight trajectories. In addition, we present computatioTi' times for the 
algorithm based on the capabilities of the TI9900 microcomputer. Conclusions 
are presented in Chapter 7. 


1.2 Previous Work on Aircraft Trajectory Optimization 

In this section, we trace the historical development of techniques for 
flight path optimization of high performance aircraft. We will first discuss 
the minimum-time problem in the vertical plane which has been under consideration 


2 


for over 35 years and then discuss the treatment of the horizontal plane problem 
and the three-dimensional problem which has occurred only recently. Finally, 
we sketch previous work in aircraft trajectory optimization using singular 
perturbation theory. 

Vertical -plane, minimum- time problem 

Before the development of supersonic aircraft, trajectory optimization was 
performed using the "quasi -steady" approximation in which the accelerations of 
the aircraft were neglected. Towards the end of World War II with the emergence 
of higher, performance aircraft, this assumption led to results which were less 
and less accurate. The high acceleration capability of the aircraft particularly 
along the flight path could no longer be neglected, Kaiser (1944) considered 
the total energy of the aircraft, expressed as energy height, to obtain a mini- 
mum-time climb path which took into account the longitudinal acceleration of the 
jet interceptor. The solution could be obtained graphically in terms of energy 
without resorting to the use of the Calculus of Variations. This graphical 
approach was later used by Lush (1951), Kelley (1952), Fuhrman (1952), Garbell 
(1953), and Lush (1956) to obtain minimum- time and minimum- fuel paths with free 
boundary conditions. Rutowski (1954) developed a graphical optimization tech- 
nique which yielded the theoretical Rutowski energy climb path. 

Much of the initial work in obtaining numerical results using the indirect 
method of the Calculus of Variations was done at RAND around 1949- 1951. Miele 
(1950, 1955a, 1955b, 1958, 1962) was active in the field of flight path optimi- 
zation for over ten years during which he introduced a new method of solving 
the minimum- time climb problem using Green's theorem. This method is applicable 
only to a restricted class of problems which can be formulated in a linear form 
in two transformed variables. Beginning in about 1954, Miele ( 1959a, b) and 
Cicala ( 1955a, b) individually and in collaboration developed the formulation of 
the fixed end-point, vertical plane problem in terms of the Bolza form of the 
Calculus of Variations. Here they considered accelerations normal to and along 
the flight path as well as control inequality constraints and state equality 
constraints. Numerical results were obtained only for simplified cases. Kel- 
ley (1959) studied zoom climbs including consideration of both normal and longi- 
tudinal accelerations and discontinuous thrust due to afterburner burnout. 

Bfyson (1966) considered the minimum-time interception of a non-maneuvering 
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target by utilizing the technique of reducing the state space by the use of 
dimensionless variables. 

In the 1960 's much emphasis was placed on the development of gradient 
algorithms as it was recognized that digital computers were necessary to solve 
flight path optimization problems. The first successful programs were developed 
by Bryson (1962) and Kelley (1962). With this technique, an optimum flight path 
is determined by comparison of an existing trajectory with its predecessors. 

The method can be made as exact as the model of the aircraft which is being 
simulated. In general, however, solutions are slow to converge, and will often 
converge on a local minimum rather than the true minimum. Balakrishnan (1969) 
proposed a modified gradient approach designed to minimize the large computation 
times. 

Development of numerical techniques for the integration of Euler-Lagrange 
equations was accomplished by Heerman (1964) and Vincent (1966). The results 
are generally 1n the form of a flooded region of trajectories where some re- 
finement is required to arrive at the desired solution. Programs of this type 
are characterized by instabilities and extreme sensitivity to particular para- 
meters . 

There was also a renewal of interest in the energy-state approximation. 

Boyd and Christie (1965) worked with the concept of energy management and de- 
veloped operational guidelines without resorting to an indirect method of the 
Calculus of Variations or gradient solution. Bryson, Desai and Hoffman (1969) 
presented a fairly complete treatment of the energy-state approximation and 
applied it to a series of vertical -plane flight path optimization problems to 
a given range. Other examples of recent applications of the energy-state ap- 
proximation to the vertical plane problem are by Meier ^ (1970), Schultz 

and Zagalsky (1972) and Parsons (1972). Parsons considered the minimum-time 
transit of a supersonic aircraft to a point which is far enough away that there 
is a central cruise arc at the maximum Mach number of the aircraft. 

Sederstrom (1972) presents results of energy management flight tests con- 
ducted with an F-8 aircraft. He showed that it was possible to calibrate an 
individual aircraft and optimize its minimum- time, energy-climb performance on 
the basis of a relatively simple procedure. Sederstrom et al . (1971) also 
consider the problem of displays and pilot workload by developing a hybrid simu- 
lation of the F-4 aircraft following optimal flight paths obtained using energy 
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management techniques. Another set of results of flight tests conducted on the 
F-8D airplane are presented by Capt. Bryan and Allison (1972). Here three 
flight. path schedules were compared; an optimum energy flight path, a flight 
manual flight path and an optimum energy flight path schedule based on nominal 
aerodynamic and performance data for the F-8D airplane. More recently. Barman 
and Erzberger (1976) considered the problem of determining optimum trajectories 
with a range constraint using the energy-state method for short-haul subsonic 
aircraft. (Jehara, Stewart and Wood (1978) investigate minimum- time loop man- 
euvers, i .e. , maneuvers in the vertical plane in which the flight path angle 
increases monotonically from 0 to 360 degrees. 

Horizontal -plane, minimum- time problem 

The minimum- fuel problem was considered in the initial work in determining 
optimum flight paths in the horizontal plane. Connor (1967) studied the sin- 
gular arc portion of the minimum- fuel path at constant altitude and Bryson and 
Lele (1969) presented the full solution of this problem. Final position con- 
straints were not included in this work. Erzberger and Lee (1971) considered 
the constant-altitude, constant-velocity minimum- time solution to a point and 
to a line. Hedrick and Bryson (1971a, b) investigated constant-altitude, variable- 
velocity, minimum-time paths to a final velocity and heading without final 
horizontal position constraints. Parsons (1972) considers constant-altitude, 
variable-velocity, minimum-time flight paths to a final point or onto a final 
line when the flight path is long enough that a cruise period at maximum velocity 
or a straight bank angle chatter arc is present in the flight. Hoffman and 
Bryson (1971) considered the case when the cruise period does not exist. 

Three-dimensional, minimum-time problem 

In 1970, Kelley and Edelbaum (1970) considered three-dimensional, minimum- 
time flight paths using the energy-state approximation and suggested an asymp- 
totic expansion procedure based on singular perturbation theory to correct the 
solution near altitude transitions. Horizontal plane final position constraints 
were not considered and numerical results were not presented. Kelley (1971a) 
further developed the asymptotic expansion approach suggested above. Numerical 
results for minimum-time paths without final position constraints were pre- 
sented by Kelley and Lefton (1972a) and Kelley (1973c). 

Hedrick and Bryson (1971, 1972) also treated the three-dimensional. 
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mini mum- time problem without horizontal-plane final position constraints during 
this same period.. Hedrick obtained more complete numerical results than Kelley 
and associates but considered less realistic flight envelope constraints. ,, 
Others have considered specific three-dimensional turns with final position 
constraints without attempting to provide a general characterization of these 
maneuvers, e.g., Stein ^ (1967), Cambell and Hartsook (1972). Parsons and 

Bryson (1972) use the energy-state approximation to consider three-dimensional, 
minimum- time flight paths to a final point or onto a final line when the flight 
path is long enough that a central cruise period at maximum Mach number is 
present in the flight. Hoffman and Bryson (1973) extend Parsons' work to con^ 
sider techniques for real-time on-line optimum flight path control using the 
reduced-order model obtained from the energy-state approximation. They also 
studied short-range maneuvers where the cruise period is absent. 

Singular Perturbation Theory (SPT) approximation of optimal aircraft trajectories 

In a series of papers appearing in the early seventies, Kelley applied the 
asymptotic expansion methods of singular perturbation theory to aircraft op- 
timization problem. By considering boundary layer correction terms he was able 
to improve the usual reduced order energy approximation in the regions' where 
instantaneous altitude transitions occur. In the first paper of the series, 
Kelley and Edelbaum (1970) considered three-dimensional maneuvers, both energy 
climbs and ^energy turns. In subsequent papers, Kelley (1970a) considered the 
general theoretical problem for a two-state system, and Kelley (1970b) applied 
the method to horizontal plane control of a rocket in a vacuum. The papers, 
Kelley (1971a) and Kelley and Lefton (1972a), consider energy state models with 
turn. More generally, Kelley (1971b, 1973c) considers three-dimensional man- 
euvers with variable mass. Note that Kelley (1973c) gives a detailed account 
of the singular perturbation approach to aircraft problems and includes most of 
the earlier work in this paper. 

More recently, other investigators have applied asymptotic techniques to 
aircraft trajectory optimization. Ardema (1976) applied the method of matched 
asymptotic expansion, one of many singular perturbation methods, to the vertical 
plane minimum time-to-cl imb problem. He calculated the zero and first order 
SPT approximations and compared them to the energy state approximation and the 
solution obtained by the method of steepest descent (which one could assume to 
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be optimal). He found that the first-order SPT approximation was close to the 
steepest descent solution and SPT required much less computation than the latter. 
In a later paper Ardema (1978) considered a general third order nonlinear SPT 
problem and studied the occurrence of singular arcs in the solution. Breakwell 
(1977j 1978) considered the vertical plane, minimum- time problem where drag D 
is much less than lift L and defined a natural perturbation parameter e=D/L. 

He also considered the occurrence of singular arcs in the solution. 

Note that the work mentioned so far only applies SPT to off-line aircraft 
trajectory optimization. In a recent series of papers Calise has applied com- 
plete time scale separation to obtain feedback controls by means of SPT. 

Aggarwal, Calise and Goldstein (1977) consider the vertical plane, minimum- fuel 
problem for a transport aircraft. Calise (1977a, 1978b) considered the vertical 
plane minimum-time problem. Calise (1977b) considered feedback control of a 
missile in the horizontal plane. Calise (1978a) considered the vertical plane 
problem to minimize a weighted combination of fuel and time for both transport 
and missile. 

This project has emphasized on-line trajectory optimization for aircraft 
control. The theoretical results (Chapter 4) address the problem of applying 
SPT to obtain feedback controls which can be computed on-line and stored on- 
board. We justify and extend the method of complete time scale separation of 
Calise (1978b) and indicate when the algorithm yields a well-defined control 
law. (See Subsection 4.3.2 and Appendix 4.2.) We also indicate methods for 
applying SPT approximations when there is no complete time scale seapration 
(see Subsections 4.3.3 and 4.3.4), i.e., the use of suboptimal solutions of the 
slow reduced order and linearization of the fast subproblem around the reduced 
order solution. Note that when the linearized fast subproblem itself exhibits 
time scale separation, one can apply the Generalized Multiple Time Scales (GMS) 
methods of Ramnath and Sandri (1969) to obtain further computational efficiency. 
Ramnath and Sinhu (1975) applied this method to determine space shuttle re-entry 
paths requiring minimum mass of heat shielding. In addition, we consider the 
state space dependence of the accuracy of the SPT approximation and the break- 
down of SPT near the terminal target (see Section 4.4). 

The practical results (Chapters 5 and 6) apply SPT to obtain a real-time 
feedback law for the three-dimensional minimum time-to-interception problem for 
a realistic aircraft model. The six-state point mass model uses real data for 
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inn II 


the aerodynamic coefficients and realistic controls and constraints (see 
Chapter 2). In addition, we obtain the SPT algorithm for this model and assess 
its real-time capability on a TI9900 microprocessor (see Chapter 6). 
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CHAPTER 2 : . 

AIRCRAFT MODEL ‘ 

2.1 Introduction 

In this chapter we describe the aircraft model used in this project. In 
summary, this model assumes a point mass approximation— that is, state variables 
describing the vehicle attitude are either omitted or used as control variables. 
This assumption implies, of course, that the characteristic time constants of 
the control system are significantly less than the time constants of the motion; 
that is, the attitude necessary for generating a certain command force or moment 
is effectively achieved instantaneously. Another assumption is that the earth's 
surface is flat, and provides the initial reference system. Next, it is assumed 
that the vehicle thrust vector is always parallel to the zero lift direction. 
According to Parsons (1972) this assumption is not a serious restriction. Fin- 
ally, vehicle mass is assumed constant. The resulting six-state equations of 
motion are presented in Section 2.2. This section also presents the state and 
control constraints for the model. 

The aircraft whose characteristics were used in this project was an early 
version of the F-4 used by Bryson, Desai and Hoffman (1969) and as Airplane 1 
by Bryson and Parsons (1971). Section 2.3 describes the numerical treatment of 
the aerodynamic coefficients characterizing this plane. In addition, this sec- 
tion describes the atmospheric model used and the computer simulation of the 
aircraft dynamics. 


2.2 Aircraft Equations of Motion and Constraints 

The system of equations for the aircraft presented below is typical of the 
point mass approximation models encountered throughout the literature. There- 
fore, the derivation shall not be repeated here. The interested reader may 
refer to Parsons (1972) from whose work we have chosen our aircraft model. 

The point mass model together with our other assumptions mentioned in the 
introductory Section 2.1 result in a six-dimensional system of first order 
nonlinear differential equations: 
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(2.1) 

X = V cosBcosy 


(2.2) 

y = V sinBcosy V TvT' .T f ; 


(2.3) 

fi = V siny ■ >0 


(2.4) 

^ _ T cosa ~ D - mg siny 
^ m 


(2.5) 

• (L + T sinot)sina 

^ mV cosy 

rr- 

.> 




(2.6) 

• ' (L + T sind)cosa - mq cosy 

^ ■ mV 



In equations (2.1) ■> (2.6), the state variables are x, y, h, V, y which rep- 
resent respectively the horizontal position (x,y) in the flat earth's plane, 
the height h above a fixed ground height, the magnitude V. of the velocity, the 
heading angle 3 in the horizontal plane of the earth and the flight path angle y* 
The reader should refer to Figure 2.1 to see the geometric relation between the 
state variables. The parameters T, D and L are respectively the thrust, drag 
and lift forces on the plane. In terms of aerodynamic coefficients these are 
given as follows: T ■ 5 - 


( 2 - 7 ) T=uT^^^(M.h) 

where T^^^ is the maximum thrust for height h an;d Mach number M. Mach number 

maX 

is related to the velocity V by, the equation 

(2.8) V = c(h)M 

where c(h) is the speed of sound at altitude ,h. 

(2.9) 

(2.10) L = L^a ’ 

(2.11) L = C|_ qS 

a 

1,2,12) Dq = Cp^qS 

(2.13) q=%pV^ 

10 ' 
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Figure 2.1 

Geometric Relationship of State and Control Variables 


In the above equations (2.9) - (2.13), is the lift coefficient, is the 

zero lift drag coefficient and n is the. aerodynamic efficiency factor. The 
parameters C|^ , Cp and ri each depend only on Mach number M. The parameter S 

denotes a characteristic surface area (49.239 m^ for this aircraft) and p(h) 
is the density at altitude h. 

The parameters u, a and o are considered control variables. The control u 
represents throttle value and varies between 0 and 1, a denotes angle of attack 
and a denotes the bank angle of the plane. See Figure 2.1 for an illustration 
of the geometric relationship of a and a to the state variables. 

In addition to the dynamic equations (2-1)- (2.6), we have the state and 
control constraints 


(2.14) 

’0 

< a 

< 

“s 

(2.15) 

0 

< u 

< 

1 

(2.16) 

0 

< V 

< 

max ' 

(2.17) 

-a 

max 

< 

0^0 

max 

(2.18) 

- 

min 

£ h 

1 


where is the stall value for the angle of attack (assumed 12° i.n our problem) 

and a is the maximum bank angle (assumed about 76° to correspond to a maxi- 

Iliu A 

mum normal load of 4g during horizontal turns). These constraints represent 
structural and controllabil ity limitations on the aircraft. 


2. 3 Numerical Simulation of Aircraft 

As we mentioned before in Section 2.1, the aircraft whose characteristics 
were used in this project is an early version of the F-4 as used by Bryson, 

Desai and Hoffman (1969) and as Airplane 1 by Bryson and Parsons (1971). The 
aerodynamic coefficients for this aircraft are tabulated as a function of Mach 
number and the maximum thrust as a function of Mach number and altitude in 

5 

Appendix 2.1. The weight of the aircraft was taken as constant at 1.5569 x lo N 

2 

and the aerodynamic reference area S is 49.239 m . The two atmospheric variables 
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required in the simulation are sonic speed and air density. , These are indicated 
in Appendix 2.1 as functions of altitude. 

The input data describing the atmosphere and the example aircraft as pre- 
sented above is in tabulated form. A continuous numerical representation of 
this data is essential for the simulation and the associated partial derivatives 
are required for the solution of the optimization problem. Therefore, the data 
was modeled using cubic spline fits which provided continuous values of the 
data with continuous first and second derivatives. These spline fits are de- 
scribed in more detail as follows: 

(a) Sonic Speed c. The sonic speed as presented in Table 3 in Appendix 2,1 

4 

is constant above an altitude of 1.2192 x 10 m. Hence a cubic spline fit was 
constructed for values of height less than the above with an end-condition of 
zero first derivative at this height. For a value of the first derivative at 
0 m, a natural* spline fit was constructed and the value of this derivative 
obtained was used as the end-condition at 0 m to construct the chosen cubic 
spline fit. 

(b) Atmospheric Density p . Since no special end-conditions were required 
to model atmospheric density, a natural spline fit was constructed. 

(c) Aerodynamic Coefficients C |^ , , n . As in Table 1 of Appendix 2.1, 

the aerodynamic coefficients are constant for M< 0.8 which gives end-conditions 
at M = 0.8 of zero first derivatives. The problem of end-conditions at M=2.0 
was solved by constructing natural spline fits which provided values for first 
derivatives at M = 2.0 that were used to obtain the required spline fits. 

(d) Maximum Thrust Table 2 of Appendix 2.1 shows maximum thrust as 

a function of two variables: altitude and Mach number. To numerically model 

this data so as to have available continuous values with partial derivatives, 
a natural bicubic spline fit was constructed. The natural fit provides zero 
second partial derivative with respect to Mach number along the altitude boun- 
dary, and at the corners zero mixed second partial derivatives. Some inter- 
polated values are presented in Table 4 of Appendix 2.1. Note that the maximum 
velocity constraint does not allow maximum thrust values at low altitude and 
high Mach number. 

The dynamic simulation of the aircraft was carried out by numerically 

*A natural spline fit provides continuous first derivatives and zero second 
derivatives at the end-points. 
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integrating the differential equations (2.1) -(2.6) using a second-qrder Adams- 
Bashforth integration routine. 
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APPENDIX 2.1 

CHARACTERISTICS OF EXAMPLE AIRCRAFT AND ATMOSPHERE MODEL 


Table 1 

Aerodynamic Coefficients for F-4 Model as Functions of Mach Number 


M 


Cl 

a 

n 

< 0.8 

0.0130 

3.44 

0.540 

0.9 

0.0140 

3.58 

0.750 

1.0 

0.0310 

4.44 

0.800 

1.1 

0.0388 

3.88 

0.830 

1.2 

0.0410 

3.44 

0.850 

1.3 

0.0408 

3.20 

0.875 

1.4 

0.0390 

3.01 

0.890 

1.5 

0.0372 

2.84 

0,910 

1.6 

0.0360 

2.68 

0.920 

1.7 

0.0354 

2.55 

0.930 

1.8 

0.0350 

2.44 

0.940 

1.9 

0.0348 

2.34 

0.945 

2.0 

0.0346 

2.25 

0.950 


-15 




0 

1.524 

3.048 

4.572 

6.096 

7.620 

9.144 

10.668 

12.192 

13.716 

15.240 

16.764 

18.288 

19.812 

21.336 


12.60* 

11.21 

9.74 

8.32 

7.07 

5.96 
4.98 
4.05 , 
3.25 

2.54 

1.96 
1.42 
0.98 
0.58 ' 
0.31 


13,70 

12.10 

10.59 

9.12 

7.70 

6.54 
5.47 

4.54 
3.60 
2.89 
2.18 
1.65 
1.11 
0.71 
0.40 


15.35 
13.48 
1,1.83 

10.36 
8.81 
7.47 
6.27 
5.20 
4.18 
3.38 
2.49 
1.87 
1.33 
0.85 
0.49 


16.86 

15.26 

13.52 

11.92 

,10.36 

8.81 

7.47 

6.18 

4.98 

3.96 

3.02 

2.27 

1.65 

1.07 

0.62 


16.06 

16.90 

15.52 

13.92 

12.14 

10.50 

8.94 

7.34 

5.96 

4.76 

3.69 

2.85 

2.05 

1.33 

0.76 


16.28 

17.13 

16.06 

14.06 

12.50 

10.77 

8.81 

7.21 

5.74 

4.45 

3.43 

2.54 

1.69 

0.98 


12.50 

10.40 

8.59 

6.81 

5.29 
4.09 
3.02 
2.05 

1.29 


11.70 

9.65 

7.70 

5.92 

4.54 

3.38 
2.31 

1.38 


*al] values times 10^ N 


Altitude 

Sonic Speed [ 

10^ m 


0 

340.2 

1.524 

334.4 

3.048 

328.3 

4.572 

322.2 

6.096 

316.1 

7.620 

309.7 

9.144 , 

303.2 

10.668 

296.6 

L 


Table 3 

Sonic Speed and Density as Functions of Altitude 


Air Density 


Kgm/m^ 


0.9045 

0.7710 

0.6525 

0.5489 

0.4583 

0.3968 



A1 ti tude 

Sonic Speed 

Air Density 

10^ m 

m/sec 

Kgm/m^ 

12.192 

295.1 

-0.3015 

13.716 

295.1 

0.2371 

15.240 

295.1 

0.1865 

16.764 

295.1 

0.1466 

18.288 

295.1 

0.1153 

21.336 

295.1 

0.07133 

24.384 

295.1 

0.04410 
























Table 4 

Maximum Thrust Interpolated Using Bicubic Splines 


xl 0 %\ 

.4 

.5 

.6 

.7 

.8 

B 


1.1 

B 



1.5 

m 


m 


m 

0 

12 . 60 * 

13.12 

13.70 

14.46 

15.35 

16.32 

16.86 

16.64 

16.06 

_ 

_ 






- 

0.762 

11.92 

12.37 

12.90 

13.57 

14.41 

15.26 

16.06 

16.64 

16.68 

- 

" 

- ■ 

- 

-■ ■ ■ 

■ - ■ 

- 

■- 

1.524 

11.21 

11.61 

12.10 

12.72 

13.48 

14.28 

15.26 

16.32 

16.90 

16.64 

16.28 

- 

■ - 

. - 

- 

- 

- 

2.236 

10.50 


11.34 

11.92 

12.63 

13.43 

14.41 

15.52 

16.37 

16.64 

16.86 

- 

- 

- ■ 

- ■ 

- 


3.048 

9.74 

10.14 

10.59 

11.17 

11.83 

12.63 

13.52 

14.55 


16.37 

17.13 

- 

- 

- 

- 

- 

- 1 

3.810 


9.39 

9.86 

10.41 

11.12 

11.88 

12.72 

13.66 

14.72 

15 . 84 . 

16.81 

- 

- 

- 


- 

- 

4.572 

8.32 

8.67 

9.12 

9.70 

10.36 

11.12 

11 . 92 - 

12.86 

13.92 

15.08 


16.68 

17.21 

- 


- 

■ ■- 

5.334 , 

7.70 


8.36 

8.94 

9.56 

10.32 

11.16 

12.05 

13.03 

14.10 


15.84 

16.55 



- 

- 

6.096 


7.34 

7.70 

8.18 


9.56 

10.36 

11.25 

12.14 

13.08 

lEia 

14.99 

15.88 

■ - 

- 

• - 

- ■ 

6.858 

6.49 

6.76 


7.56 

Iwnl 

8.81 

9.56 

10.41 


12.23 

13.26 

14.23 

15.08 

- 

- 

- 

- 

7.620 

5.96 

6.23 

6.54 

6.94 

7.47 

8.10 

8.81 

9.61 

10.50 

11.48 

12.50 

IHEO 

14.23 

14.81 

15.39 

- 


8.382 

5.47 

5.69 

6.01 

6.36 

6.85 

7.43 

8.14 


9.74 

10.68 

11.70 

IBB 

13.39 

14.01 

14.63 

- 

■- 

9.144 

4.98 

5.20 

5.47 

5.83 

6.27 

6.85 

7 . 47 - 

8.18 


■ 9.83 

10.77 

11.65 

12.50 

13.21 

13.83 

- 

- 

9.906 

4.49 

4.72 

4.98 

■ 5.34 

5.74 

6.23 

6.81 

7.43 

8.14 

8.90 

9.79 

10.63 

11.48 

12.19 

12.81 

- 

- 

10.668 

4.05 

4.27 

4.54 

4.85 

5.20 

5.65 

6.18 

6.72 

7.34 

8.05 

8.81 

9.61 

10.40 

li . l 2 


12.19 

12.63 

11.430 ■ 

3.65 


4.05 

4.31 

4.67 

5.07 

5.56 


6.63 

7.25 

7.96 

8.72 

9.47 


EBB 

11.12 

1 !^ 

12.-192 • 

3.25 

3.43 

3.60 

3.87 

4.18 

4.54 

4.98 


5.96 

6.54 

7.21 

7.92 

8.59 

9.16 • 

9.65 

MiMi a 

10.45 

12.954 

2.89 


3.25 

3.47 

3.78 

4.09 

4.45 


5.34 

■ 5.87 


7 . 07 - 

7.70 

, 8.23 

8.67 

Kff 2 

9.39 

13.716 

2.54 

2.71 

2.89 

3.11 

3.38 

3.65 

3.96 

4.31 

4.76 

5.25 

5.74 

6 . 27 ‘ 

6.81 

7.30 

7.70 

8.05 

8. 36 

14.478 

2.22 

2.36 

2.54 

2.71 

2.94 

3.20 

3.47 

3.83 

4.18 

4.63 


5.56 

6.01 

6 . 41 - 

■ 6.76 

7.07 

7.34 

wmm 

1.96 

2.05 

2.18 

2.31 

2.49 

2.71 

3.02 

3.34 

3.69 

4.05 

4.45 

4.89 

5.29 

5.65 

5.92 

6.18 

6.41 


1.69 

1.78 

1.91 

2.00 

2.14 

2.38 

2.62 

2.94 

3.25 

3.56 

3.91 

4.31 

4.67 

4.94 

5.20 

5.38 

5.56 

16.764 

1.42 

1.56 

1.65 

1.73 

1.87 

2.05 

2.27 

2.54 

2.85 

3.11 


3.78 


4.36 

4.54 

4.72 

4.85 


1.20 

1.29 

1.38 

1.47 


1.73 

1.96 

2.18 

2.45 

2.71 

2.93 

3.29 

3.56 

3.78 

3.96 

4.09 

4.23 

13.288 

0.98 

1.02 

1.11 

1.20 

1.33 

1.47 


1.82 

2.05 

2.27 

2.54 

2,80 


3 . 20 , 

3.38 

3.51 

3.60 

19.050 

0.76 


0.89 

0.98 

1.07 

1.20 

1.33 

1.51 

1.69 

1.87 


2.31 

2.54 

2.71 ‘ 


2.94 

2.98 

19.812 

0.58 

0.67 

0.71 

0.76 

0.85 


1.07 

1.20 

1.33 

1.51 

1.69 

1.87 

2.05 

2.18 : 

2.31 

2.36 

2.40 

wsaaim 

0 . 4,5 

0.49 

0.53 

0.58 

0.67 

0.76 

0.85 

0.93 

1.02 

1.16 

1.33 

1.51 

1.65 

. 1.78 

1.82 

1.87 

1.87 

21.336 

0.31 

0.36 

0.40 

0.44 

0.49 

0.53 


0.67 

0.76 

0.85 


1.16 

1.29 

1.38 

1.38 

1.38 

1.42 


*all values times 10** N 


































CHAPTER 3 

THE OPTIMAL CONTROL PROBLEM AND EXACT, ITERATIVE METHODS OF SOLUTION 
'Ll Introduction 

In this chapter we describe the trajectory optimization problem considered 
in this project and some exact, iterative methods for its solution. We wish 
to draw particular attention to the continuation methods described in Subsection 
L3.3. These methods seem to offer very efficient trajectory optimization al- 
^or ithms, although they are unfortunately not yet fast enough for on-line optimal 
control of aircraft. 

The optimization problem chosen for this problem was the minimum time in- 
terception problem. That is, the problem was to minimize the total trajectory 
time from a given initial point in the six-dimensional state space to a given 
final point. In Section 3.2 we describe the mathematical formulation of this 
optimization problem. 


;■ . 2 F ormulation of the Minimum Time Interception Problem 

Here we formulate the minimum time interception problem as in Bryson and 
•; c (1375) or in Athans and Falb (1966). First, we have a dynamic system 

(2,1) = f(x,u) 

.'here x denotes the state vector and u denotes the control variable. We have 
underlined the state x and control u to distinguish them from the horizontal 
position x and the throttle control u. In later sections we will omit the 
underlines when the context makes clear which case is meant. In our case 
(>s,yLhsV,3,Y) and u=(u,a,a). Equation (2.1) is the vector representation 
-T'f tr.e system (2,2. 1) - (2.2.6). 

In addition to the dynamic equations (2.1) there are inequality constraints 
; 2) C(x,u) <0 

'’h;ch must be satisfied at each time of the trajectory. Equation (2.2) is the 



vector representation (hence 0 rather than 0) of the state and control con- 
straints given in (2.2. 14) - (2. 2. 18) . 

At the initial time tQ an initial state Xg is specified, 

(2.3) X(tg) = 2Sp 

where * Similarly, at the final time t^ a final state 

x^ is specified,, 

(2.4) x(t,p) = x^ 

where x^ = (x,p,y^,h^,V^,3f sYf) . Note that the initial time is fixed but the final 
time need not be. 

Equations (2.1) -(2.4) denote constraints on the possible trajectories of 
X, u. The optimization problem is to minimize the cost of the trajectory where 
the cost has the form 


(2.5) 


J(u) 



x,u) dt 


For the minimum time problem L is simply the constant. 


(2.6) L(x,u) = 1 


Let us say that a control trajectory u is feasible if u together with its 
corresponding state trajectory x, found from integrating (2.1) with initial 
conditions (2.3), satisfies the constraints (2.2) and the terminal condition 
(2.4). The optimal control problem is to find a feasible control trajectory 
u* such that 

(2.7) J(u*) < J(u) 

for all other feasible controls u. Solutions of the optimization problem are 
.usually found by solving the Euler-Lagrange first order necessary conditions. 
These conditions are expressed in terms of a Hamiltonian function defined 

(2.8) H(x,X,u) = L(x,u) + X^f(x,u) 
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where A denotes a vector of the same dimension as x and A^ denotes its transpose. 
In our case, H is given by 

(2.9) H = 1 + A V cos3cosy + A V sinBcosy 

X y 


T cosg - D - mg siny 
m 


+ A^V siny + A^ 


The Euler-Lagrange equations are a system of differential equations in x and A, 
namely 


» 

(L + T sina)sina 

+ \ 

(L + T sina)cosa - mq cosy 

mV cosy 

’ A 

y 

mV 


(2.10) ^=f(x,u) 

(2.11) ^=G(x,A) 


For optimal u*, x* there is an optimal A* such that (2.10), (2.11) are satisfied 
and such that 

(2.12) H(x*,u*,A*) = min H(x*,u,A*) 

u 

at all times. This relationship is known as the minimum principle (see Athans 
and Falb (1966)). 

In principle we can solve (2.12) for u in terms of x, A and eliminate u 
from (2.10), (2.11). The result is a system of differential equations in x and 
A only, 

(2.13) ^=F(2,i) 

(2.14) ^=G(x,X) 

with mixed initial and boundary conditions on x from (2.3) and (2.4) rather than 
initial conditions on both x and A. Thus, (2.13), (2.14) with boundary con- 
ditions (2.3), (2.4) is a two-point boundary value problem (TPBVP). Many solu- 
tion techniques solve the original optimization problem by solving this TPBVP, 
Note that when t.^ is not specified, we obtain the extra condition 
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(2.15) H(x*.X*.u*) = 0 


at all times. This extra condition can be used to find t^. 


3.3 Exact, Iterative Solutions of the TPBVP 

Whether the full system solution or the solution to one of the reduced- 
order formulations of the problem is being sought, a numerical solution to the 
TPBVP is necessary. Order reduction is usually achieved through the application 
of Singular Perturbation Theory (SPT). For a system which is considered 
"stiff"— that is, the time scales of some of its state variables are signifi- 
cantly faster than those of the remainder — the approximation arising from a 
reduction in order can be made to be a reasonably good one, if done judiciously. 

Assuming, then, that the system has been reduced to a TPBVP, the solution 
procedures considered in this project fall into three general categories: 

(i) Steepest Descent (Gradient) Methods 

(ii) Quasi linearization Techniques 

(iii) Continuation Methods Using a Parameter 

These methods will be discussed in varying detail in the following sections. 
3.3.1 Steepest Descent Techniques 

These are the most widely used methods, and their strengths and weaknesses 
are well known. Inequality constraints are typically handled by using penalty 
functions. There are often convergence problems, due to the presence of state 
variable inequality constraints and singular arcs. This is because of the ab- 
sence of control variables in the inequality constraints or in the gradient of 
the cost function. It is expected that some of these numerical problems may be 
alleviated by using a "generalized gradient method" as described in Mehra and 
Davis (1972). Briefly, this method uses the constraints to dictate, at each 
step, which of the entire set of control and state variables are to be selected 
as control variables for the next step. The gradient of the cost function with 
respect to the independent variables, called the generalized gradient, is then 
computed by solving a set of equations similar to the Euler-Lagrange equations. 
Directions of search are found using gradient projection and the conjugate 
gradient method. The procedure, then, is based on the idea that there is no 
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real mathematical distinction in the use of u(t) or elements of u(t) and x(t) 
as the independent or manipulative variables of the system x=f(x,u). Criteria 
for selecting the independent set include: any variable which lies on the con- 

straint boundary should be included in the set; choosing some of the state 
variables as the independent variables often improves the rate of convergence; 
the independent set must be chosen so as to retain recursiveness to avoid in- 
verting large matrices. 

The details of this technique are in Mehra and Davis (1972). It is an- 
ticipated that algorithms based on steepest descent, even when applied to the 
reduced TPBVP's, may be too slow for real time solution of any of the basic 
problems. This method will be useful, however, in developing the full -system 
solution. 

3.3.2 Quasil inearization Methods 

Basically, these methods revolve around doing a linearization around a 
zero-order solution. Such a solution would arise from setting the "approxima- 
tion parameter," e, to zero, as is done in SPT; or, as is done in the case of 
continuation methods (next section), e is set to some e^, possibly zero, for 
which the solution is readily obtained. The appeal of quasilinearization tech- 
niques lies in the fact that they are far less sensitive to changes in the ini- 
tial conditions than the shooting continuation methods described in the next 
section. Thus, in advancing the parameter e from Cq or 0 to its "real" value, 
generally 1, some combination of techniques based on quasilinearization method- 
ology and continuation theory may be developed. This will hopefully allow the 
exploitation of the advantages of both schemes, utilizing one where the other 
is weak. One anticipated disadvantage of the quasil inearization-continuation 
approach is the computational effort required to step the parameter e along 
to its final value. 

Basic references on the theory of quasil inearization are Bellman and 
Kalaba (1965), Dyer and McReynolds (1970), Polak (1971), and Keller (1968). As 
in Subsection 3.3.1, it is anticipated that inequality constraints can be 
handled by means of penalty or barrier functions. This should be true both for 
this section and the following one. 

Quasil inearization techniques are cast in the form of Newton-Raphson 
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problems in Polak (1971). Because of the comparison of this form to the con- 
tinuation methods of the next section, Polak *s algorithms will be presented 
here. The outline in Dyer and McReynolds (1970) is presented more directly in 
terms of the calculus of variations problem in optimal control. 

For quasi linearization, g is augmented by a third vector element. Suppose 
we have a differential equation 

^ ~ "^(x, t) 


with initial and final boundary conditions 

(3.2) = 0 

(3.3) g^(x(t^)) = 0 

Note that in the control problem, x in (3.1) would actually include both the 
state X and the adjoint X. That is, x= (x,X) and (3.1) is given by (2.13), 
(2.14). Let g(xQ,x(*)) be the function defined by 


(3.4) 


g(xQ,x(-))(t) = Xq 



and define g from g^, g^ and g so that 


(3.5) g(xQ,x(*)) 


g(xQ,x(-)) 

9o(x(to)) 

g^(x(t^)) 


Thus, g maps the pair (xp,x(*)) into the pair (yg,y(-)) where yQ = (gQ(x(tQ) ) , 
9f(x(t^))) and y(t) = g(xQ,x(* )) (t) . In the following, let L denote the linear 
space of pairs (xp,x(*)) where Xg is a vector and x(-) is a piecewise contin- 
uously differentiable function of t. Then g is a continuously differentiable 
map of L into itself. We will assume that [9g(z)/8z]”^ exists for all z in a 
sufficiently large subset of L, where z- (xg,x(*)). With this assumption, the 
Newton-Raphson method for solving g(z) = 0 is defined by 
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(3.6) 


9z ^ 


z) 


-g(z) 


where z is the iteration for the solution. Thus, "^z is found from i by 
utilizing the above relationship, exactly as in the finite-dimensional case 
(next section). Due to the presence of x(t) £ t^jt^] in z, the Newton-Raphson 
method above has been formulated in a Banach space L. 

Substituting the g into the above equation, rearranging terms, and dif- 
ferentiating the terms associated with g with respect to time, there results 


(3.7) 


^ 3 . f . (x(t).t) + f(x{t),t), 

j+1 j+1 

9o( X (tg)) = 0, g^( x (tf)) = 0 


t e[tQ,tf] 


This differential system is called the quasi linearization version of the Newton- 
Raphson method. McGill and Kenneth (1963) and Bellman and Kalaba (1965) provide 
more details. The algorithm proceeds as follows: 

1) Select an x(-) £ C^f^jCtQjt^] such that gQ = g^=0; if Xq and x^ are such 
that the boundary conditions are met, then an acceptable x(t) may be 


x(t) = Xp + [(t-tQ)/(t^-tQ)](x^-XQ), te[tQ,tf]; 

2) Set j = 0; 

3) For te [tQ,t^], compute f(x(t),t) and ~ (x(t),t); 

4) Compute x (t) by solving the Newton-Raphson differential system, in- 
tegrating in a stable direction (or else, combining the technique of 

Roberts and Shipman (1967), described in the next section); 
j+1 j 

5) If X (t) is "close to" x(t) by some standard, stop; 

6) Otherwise, set j = j+1 and go to step (3). 

Quasi 1 inearizati on techniques such as the one outlined above provide at 
each step an approximation to a solution of x=f(x(t),t) which satisfies the 
boundary conditions. This is a major difference between quasilinearization 
techniques and procedures outlined below. The latter, at each step, provide 
an approximation to the solution and the boundary conditions which satsify the 
differential equation. 
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3.3.3 C^Q^n t i n ua t i on Methods f or the Solution of the TPBVP 

In this section we are interested in describing some continuation or im- 
bedding methods for solving the two-point boundary value problem and discussing 
their advantages over more conventional methods. To this end, we have omitted 
many mathematical details necessary for a logically rigorous presentation of 
this material. We hope this omission will make the presentation clearer for 
the reader unfamiliar with continuation methods. Mathematically rigorous results 
may be found in Ortega and Rheinboldt (1970). 

For definiteness, consider the following two-point boundary value problem 
(TPBVP) which depends on a scalar parameter: 

TPBVP : Find x(t,e), a(e) and t(c) for times t and all parameters e such 
that x(t,e), a(e) are vectors in and t(c) is a scalar time, and such that 
these quantities satisfy 

(3.8) If (t,e) = f((t,£)t,£) 

(3.9) x(0,e) = a(e) 

(3.10) (j)(a(e),x(T(e),e),T(e),E) = 0 

for all e and all t with 0 < t < t(e) . 

Note that the function (j), which maps R^xR'^xR^xR^ into is given 

beforehand, and equation (3.10) summarizes the n+1 initial and final conditions 
necessary to deduce n initial conditions (the vector a(e)) and the terminal time 
x(e). One can choose (p so that (3.8) - (3.10) represent almost any initial or 
boundary value problem. In particular, (3.8)- (3.10) can model the TPBVP for 
which the terminal time t is not given explicitly (this is the situation for 
minimum time control problems). 

To solve the TPBVP we transcribe the equations (3.8)- (3.10) to a system 
of n+1 nonlinear equations for the initial condition a(e) and the final time 
t(e). To make this transcription define the function x(t;a,e) as the solution 
of the initial value problem 

(3.11) || (t;a,e) = f(x(t;a,e) ,t,e) 

(3.12) x(0;a,e) = a 
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for all t, a, e. Then we solve for a(e) and T(e) from the equation 

(3.13) (})(a(e),x(T(e); a(e) ,e) ,T(e) ,e) =0 

Note that (3.13) represents n+1 equations which we desire to solve for the 
initial state a(e) (n conditions) and the final time x(e) (one more condition). 
The TPBVP of (3.8)- (3.10) has a solution if and only if (3.13) has a solution. 
Thus, the TPBVP is reduced to solving the nonlinear equation (3.13). 

Different techniques for solving the TPBVP derive from techniques for 
solving the equation (3.13). Thus, let us define the function G(a,x,e) as 

(3.14) G(a,x,e) = c{)(a, x(x;a,e), x, e) 

and let us generically represent (a,e) by the n+1 vector v. Then (3.13) takes 
the more general form 

(3.15) G(v,e) = 0 

which we solve for v as a function of e. 

Let us suppose that e varies between 0 and 1. Often our problem is to 
solve a difficult problem 

(3.16) G(v,l) = 0 

when we know how to solve an easier problem 

(3.17) G(v,0) = 0 

Sometimes the parameter e occurs naturally in the problem (i.e., the viscosity 
in a hydrodynamic problem), but often we introduce the parameter e artificially. 
In either case, the rationale for replacing the single equation (3.16) with a 
family of equations (3.15) is that we may be able to continue the solution at 
e = 0 in (3,17) to the solution at e=l in (3.16) more easily than computing the 
solution at e=l in (3.16) by itself. For example, a classic technique for con- 
tinuing the solution v(0) of (3.17) to a solution v(e) of (3.15) for e>0 is to 
expand the vector function v(e) in a Taylor series in e. The validity of the 
expansion requires some regularity of the function G with respect to e, but in 
some cases where G does not depend regularly on e it is often still possible to 
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use perturbation analysis to find a singular perturbation expansion of v(e) in 
terms of e. 

As a simple but concrete example, consider the following two quadratic 
equations: 

(3.18) y + y^ = e 
and 

(3.19) X + ex^ = 1 

The solution at e = 0 of 
y(e) of (3. 18) for e > 0 
series expansions 

(3.20) y(e) = e - 
or 

(3.21) y(e) = -1 - e + + . . . 

This series corresponds to the Taylor expansion in powers of e of 

(3.22) y(e) = 

and 

(2.23) y(e) = 

Equation (3.19) is not regular with respect to e as one can see by setting 
e=0; there is only one solution x(0) = 1, whereas there must be two whenever 
Nevertheless, one can still expand the solution x(c) of (3.19) in the 
singular perturbation series 

(3.24) x(e) = 1 - e + . . . 
or 

(3.>25) x(e) = -1 + 2+ ... 

e 

which series correspond to the exact solutions 


equation (3.18) gives y(0) = 0 or y(0) = -l. The solution 
depends regularly on e and has the two possible power 

+ . . . 
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(3.26) x(e) = 


and 

(3.27) x(e) = 
respectively. 

The approximations (3.20), (3.21) to (3.18) or (3.24), (3. 25) to (3. 19) are 
good when e is "small," but the trouble is that 1 is not small. In fact, 
the expansions (3.20), (3.21) and (3.24), (3.25) are valid only when \e\<h. 

For |el > %, the series do not converge and one cannot use finitely many terms 
of the series to approximate the exact answer. For example, if one tries to use 
the series from (3.20) for e=l to approximate (3.22), one obtains the approxi- 
mations y(e) - 0,1,-1,2,-5,14,-42,132 and so on, which become progressively 
worse as one adds more terms to the series in (3.20). The exact answer is (3.22) 

with 1, which gives y(e) = .618033989 

In our examples in (3.13) or (3.19) perturbation analysis cannot continue 
the e=0 solution beyond e = %, Nevertheless, there are continuation methods 
which can continue the e = 0 solution all the way to e= 1. One such method is 
the method of differentiation with respect to a parameter of Davidenko (1953). 

Consider the general equation (3.15), G(v,e) = 0, and suppose that there is 
a solution v(e) of (3.15) which is continuously differentiable with respect to 
e. Taking the derivative of equation (3.15) with respect to e gives us the 
following differential equation for v: 

(3.28) G ^ + G =0 

' ^ V de e 

In (3.28) the expression is the matrix of first partial derivatives of G with 
respect to v, and G^ is the vector of first partial derivatives with respect to 
e. The initial conditions for the differential equation (3.28) is just the v(0) 
given by equation (3.17), namely 

(3.29) G(v(0),0) = 0 

Assuming that we can solve the linear equation 
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(3.30) G^(v,e)w + 6^(v,e) =0 

for w given v and e» then we can numerically integrate equation (3.28) from the 
initial condition (3.29). If the matrix 6^(v,e) becomes singular and (3.30) 
has no solution, then the numerical integration of (3.28) may fail. However, 
Keller (1977) presents methods for continuing the integration when G^(v,e) 
becomes singular. In this case one finds that (3.28) has bifurcating solutions. 
We discuss such bifurcations further in Appendix 3.1. 

As a simple example consider equation (3.28) again. In this case (3.28) 
becomes 

(3.31) (l + 2y)^=l 

with two possible initial conditions, either y(0) = 0 or y(0) = -l. For all e^O, 
we can integrate (3.31) without difficulty and obtain the two solutions 

(3.32) y(e) = 

from the initial condition y(0) = 0, and 

(3.33) y(e) = ~ 

from the initial condition y(0) = -l. 

Note that Davidenko's method succeeds in finding the solutions at e=l 
whereas the power series method in (3.20), (3.21) fails. Using an integration 
step size of h=.l and Euler's method of integrating (3.31), we obtain the ap- 
proximation .6372 for y(l) corresponding to y(0) = 0. Compare this to the exact 
solution y(l) = .6180. . . . 

The Newton-Raphson method is an alternative for solving the equation (3.16), 
G(v,e) = 0, directly for any particular e, provided one has an initial estimate 
of v(e). Recall that the Newton-Raphson method calculates successive approxima- 
tions Vj^(e) from the recursive formula, 

(3.24) “ G^(V|^(e),e)"^G(V|^(e),e) 
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To start the method one requires an initial estimate Vj(e). In essence, if 
Vj.Ce) is "reasonably close" to v(e), then the estimates V|^(e) computed from (3.34) 
converge very quickly to v(e) as k tends to infinity. The disadvantage of using 
Newton-Raphson's method is that we may not have an initial which is reason- 
ably close to the actual solution v. A continuation method such as Davidenko's 
method can help overcome this difficulty. 

As a specific example, consider the trivial equation 

(3.35) 1 - = 0 

and suppose that we know that v(10) = 10 but we wish to approximate v(0) from 
this initial guess. The guess v^(0) = 10 is disastrous for applying Newton- 
Raphson's method. The recursion equation (3.34) for e = 0 becomes 

(3.36) = Vj^ + 1 - e^*^ 

and we find that if Vj=10, then V 2 ~ -22015, -22014, -22013 and so on. 

One requires over 22,000 iterations of (3.36) to approach the true root v = 0. 

Davidenko's method leads to the numerical integration of the differential 
equation 

(3.37) ^ = 1 

with the initial condition v(10) = 10. If we use Euler's method of integration 
with step size h, we will require roughly 10 xh"^ integration steps to reach 
the approximation of v(0). Note that the approximation will be accurate to 
order h. Thus, Davidenko's method leads to the approximate solution of (3.35) 
at e = 0 in a reasonably few integration steps if we do not require' great accuracy 
in our approximation. ' ^ 

For future reference we now present another continuation method different 
from Davidenko's method. Ortega and Rheinboldt (1970) present this method as 
an example of continuation methods. 

This continuation method differs from Davidenko's method in that it does 
not use a differential equation to calculate v(e) at successive values of e. 
Instead it uses the Newton-Raphson method (3.34) in a recursive fashion as we 
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now explain. 

Let us sol ve G(v,e) = 0 at e = e. ,e«, . . . ,e where 0 < e. < . 

± ^ ^ i\ 

Suppose that we have an initial estimate for v(0) or perhaps we know v '^ vCO).' 
Also assume given the sequence m^,m2»...»m^ of positive integers. From this 
initial estimate one defines recursively as follows. For l^k<m^ -l, 

define 

(3.38) = V|^(e^.) - G^(V|^(e,.),e^.)‘^G(V|^{e^).e^.) 


for k=l and i>l, define v^(e.j) as technique operates by 

recursively calculating Newton-Raphson approximations of the equation G(v,e..)=0 
by using as an initial estimate for v(e^) the approximation of found, 

from solving G(v,e.._j^) = 0. By suitably choosing the e.| and the m^., one can 
continue a solution of G(v,0)=0 to a solution of G(v,e) = 0 for relatively large 
values of e. Accuracy depends on how close together one takes the e- and how 
large one takes the m^. . 

The Newton-Raphson technique is a fast converging approximation method 
provided that the initial estimate is close to the exact solution. In this case, 
the error at each step is reduced by squaring. That is, the error at the 
k+1 iteration is approximately and thus, convergence is very fast. A 
continuation method such as Davidenko's method is inferior to Newton-Raphson 's 
method from the standpoint of accuracy with respect to computation speed. The 
Davidenko method makes a final error proportional to the integration step size 
h in the integration of (3.28) (or proportional to h^ for a>l if a better in- 
tegration scheme^^than Euler's method is employed). If one tries to obtain good 
accuracy by choosing h small, the integration time may be quite long. For ex- 
ample, numerically integrating (3.37) to obtain the solution v(0) of (3.35) to 
four decimal places accuracy would require about 10^ integration steps using 
Euler's method of- integration. 

However, although the continuation methods are inefficient in computing 
highly accurate solutions, these methods are much less sensitive than Newton- 
Raphson to the initial estimate of the solution. 

Similarly, continuation methods are not as efficient as perturbation ex- ^ 
pansions when the perturbation parameter is small. On the other hand, continuation 
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methods such as Davidenko's method provide solutions even when the perturbation 
parameter is large. 

Ideally, one might use a continuation method together with a method such 
as Newton-Raphson or perturbation analysis. For example, one might first use 
Davidenko's method with a modest integration step size to compute an approxima- 
tion reasonably close to the exact solution. Then, using Newton-Raphson 's 
method and using the Davidenko approximation as an initial estimate, one could 
obtain a very accurate approximation of the exact solution. Such an algorithm 
is described in Appendix 3.2. 

Having presented these different continuation techniques for computing the 
solution of G{v,e) = 0, we now show how to apply these techniques to the TPBVP. 
There are basically two continuation methods for solving the TPBVP which corres- 
pond to the two continuation methods for solving the equation G(v,e) = 0. Kubicek 
and Hlavacek (1973) use Davidenko's method to solve the TPBVP, and Roberts and 
Shipman (1967) use the continuation version of Newton-Raphson to solve TPBVP. 


Kubkcek and Hlavacek Algorithm 

They also call this method general parameter mapping or GPM. The function 
G is given by equation (3.14) and then the TPBVP is solved by solving (3.15). 
Kubicek and Hlavacek do this by using Davidenko's method--that is, by solving 
the differential equation (3.28). To use Davidenko's method, we must compute 
G and G as follows: 

V E 


(3.39) 


G^(a,T,e) 


‘•’x • "'a \ 
Vx ' \ * 


and 


(3.40) G^(a,T,e) = <p^ • + (p^ 


The subscripts in equations (3.39) and (3.40) denote partial derivatives. Thus, 
(|) =1^ is the (n+1) X n matrix of first partial derivatives of (j)(a,x,c}),e) with 
respect to a (remember that p takes its values in R" and a is a vector of R^). 

Similarly, p is the (n+l)xp matrix of derivatives of cj) with respect to x; x 

X d 

is the nxn matrix of derivatives of x(T;a,e) with respect to a. The expres- 

3 X 3 X 

sions X and x denote the n-vectors ^ and — . Likewise, p and (|) denote the 

T 0 oT O0 T 0 

(n+l)-vectors and 
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From the original formulation of the problem, one can find <P> <p , <p 
and <t>^. However, to compute and 6^, one needs x, x^, x^- and x^., and one 
must integrate (3.11), (3.12) directly. That is, one integrates (3.11), (3.12) 
to t = T to obtain x(T;a,e). Then x^(x;a,e) is given from (3.11) as 

(3.41) x^ = f(x(T;a,e),T,e) 

To obtain x^ and x we differentiate (3.11), (3.12) with respect to a and e re- 
a c 

spectively. Thus, x^ is the solution of the initial value problem 
9x 

(3.42) -^ = fj^(x(t;a,e),t,e)-x^ 

(3.43) x^(0;a,e) = I 

where I is the n^n identity matrix. Likewise, x is the solution of the initial 
value problem 

dx 

(3.44) 3 ^ = f^(x(t;a,e),t,e) *x^ + f^(x(t;a,e) , t,e) 

(3.45) x^(0;a,e) = 0 

Thus, an algorithm for solving the TPBVP (3.8), (3.9), (3.10) for 0<e<l 
might be the following: 

1. Solve the problem for e = 0. Set e<0 and a(e)=aj^, t(£) = Tq. Choose 
integration step size h (for integration with respect to e). 

2. For e, a(e), x(e) given, compute x, x^, x^ and x^ from integrating 
(3.11), (3.12), (3.41), (3.42) - (3.45). 

3. With these values of e, a(s), x(e) and x, x,, x , x from #2, use 

a X G 

(3.39) and (3.40) to compute and G^. 

4. Using some numerical integration scheme on (3.28), compute 

^ and use this to evaluate (a(e+h), x(e+h)). 

5. Update e to e+h. If g>1, stop. If e<l, return to step 2. 

Roberts and Shipman 

Roberts and Shipman (1967, 1968) solve a TPBVP 
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(3.46) ^=f(x.t), Oststj 

(3.47) x(0) = a 

(3.48) x(a,x(t^)) = 0 

where is given by allowing the terminal time to vary from 0 to t^ and keeping 
all other conditions the same. In terms of the general formulation (3.8), (3.9), 

(3.10) , Roberts and Shipman solve the problem 

(3.49) II = f(x(t;s),t) 

(3.50) x(0,e) = a 

(3.51) x(a,x(T(e),e) ) = 0, t(e) = e 

for all e and all t such that 0<t<T(e) = e. Thus, Roberts and Shipman allow 
just the final time t to be the parameter e, and they solve the original problem 
(3.46) - (3.48) by varying e from 0 to t^ in (3.49) - (3.51). The continuation 
method they use is the extension of Newton-Raphson' s method we described for 
solving G(v,e) = 0. Since this continuation method is easily applied to the more 
general problem (3.8) - (3. 10) , we do this rather than treat the specific problem 
Roberts and Shipman (1967, 1968) use in their papers. Nevertheless, for future 
reference we will refer to this continuation solution of the general TPBVP as 
the Roberts-Shipman method. 

The method consists of applying the equation (3.38) to solve G(v,e) = 0 when 
G is defined by (3.14), One obtains G^ from (3.39) just as in using Kubicek 
and Hlavacek's GPM method. This requires numerical integration of the equations 

(3.11) , (3.12) and (3.41) - (3.43) just as before. One possible algorithm to 
solve the TPBVP using the Roberts-Shipman method is then the following: 

1. Solve the problem for e = 0 to find ag, Tq. Choose £q = 0, 

0<e^<C2 ••• choose positive integers m^. ; l<i<n. Set 

i = 0, k = 0, a^{e.) = '^k^^i^'^O* 

2. For e.j, a^{e.) > T|^(e.j) given, compute x, x^, x^ from (3.11), (3.12), 
(3.41) - (3.43). 

3. With these values of , 3|^(c.|), x, x^, x^ from #2, use 

(3.39) to compute G^. 
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4. 


Invert from #3 and evaluate 
equation (3.38) (the Newton-Raphson step). 

5. Update k to k+1. If k = m^. , go on to step 6; otherwise return to 
step 2. 

6. If i = n, stop. Otherwise set 3 q(^*}+i) ~ ^k^^i ^ ^O^^i+1^ ~ ^k^^i ^ 

and reset k to 0; set i to i+1 and return to step 2. 

Summa ry of H isto ry and Advantages of Continuation Methods 

Continuation or imbedding methods have long been used to prove existence 
theorems for operator equations. Ortega and Rheinboldt (1970) give a nice 
discussion of the method with many historical notes and references. Ficken (1951) 
contains references and notes on the literature previous to 1950, including some 
from the last century. 

Lahaye (1934) and later, independently, Davidenko (1953) first applied con- 
tinuation methods to the numerical solution of nonlinear equations of the form 
6(v,e) = 0. Davidenko introduced the method of differentiating G(v,e) with re- 
spect to the parameter e to obtain a differential equation for v. 

More recently, the method of continuation has been used to solve numerically 
fixed point problems. See Kellogg, Li and Yorke (1976) for example. Rigorous 
and powerful mathematical treatments of the continuation method rely on topolo- 
gical homotopy theory. 

A preliminary study of simple nonlinear problems indicates the method of 
differentiation with respect to a parameter, Davidenko 's method, when used in 
conjunction with perturbation analysis or a Newton-Raphson method, offers a 
powerful numerical technique for solving nonlinear problems. The main advan- 
tage of the continuation method is that it permits one to approach an exact 
solution* even when the initial estimate is not close to an exact solution. If 
one uses a method such as Newton-Raphson together with the continuation method, 
then one can also achieve good computational efficiency. The essence of the 
technique is to use the continuation method to obtain a rough first estimate 
and then use Newton-Raphson to refine the estimate. See Rheinboldt (1978). 

Many other variations are possible. For example, Keller (1976, “Bifurca- 
tion Theory and Nonlinear Eigenvalue Problems," unpublished lecture notes, Cal- 
tech) poses a continuation problem in terms of an equation G(v,e)=0 where v is 


*This includes all local minima which satisfy the first order necessary con- 
ditions of optimality. 
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an element of a Banach space and G(»,e) is an operator on that Banach space. 

Such an approach allows one to treat the TPBVP as a nonlinear operator equa- 
tion on the Banach space of solution functions x(»,e) of (3.8) - (3.10). Al- 
though seemingly more complicated than the approach of Kubicek-Hlavacek or 
Roberts-Shipman, this infinite dimensional point of view may offer the advan- 
tage of a numerically more stable solution algorithm. 

The continuation methods seem very promising for application to the non- 
linear TPBVP which appear in optimal control problems, and the technique 
deserves wider circulation. 

Example 

To see how continuation methods might be applied to aircraft trajectory 
optimization, consider the planar minimum-time-to-cl imb problem of Ardema (1976). 
Because of the absence of terminal constraints on horizontal position, the 
variables x and y are unnecessary. Also, the problem is solved in the vertical 
plane, so that the heading, B, and roll control variables are omitted. Thus, 
h, Y V remain as state variables and velocity V is replaced by energy e. 

Also, thrust T is assumed constant, so that L (or a) is the only control variable. 

This reduced system is of a dimension six which would allow a relatively 
inexpensive means of verifying and developing the various computational al- 
gorithms. The TPBVP for this system looks like 


where 

and 


where 


X = f(x(t),n»e) 

X ~ (h,y,e,A|^,X^,A^) 

h = V siny 
Vy = L* - cosy 
e = V(T- D) 

= Xj^(g/V)sinY - X^(l/V^)(g/V)(L* - cosy) - X^ 

X - -X. V cos - X siny/V 
y h Y 

X^{l/v2)(l/V)(L*-cosy) - |£ 
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I 


p = V(T-D) 

D = D(L*^,h,e) 

eX 

L* = ^ , the optimal control. 


In the expression L*, B = B(M) is an induced drag parameter, suitably dimension- 
al i zed for the available tabular data, and e is the parameter used as the ‘ 
independent variable for solving the vector of unknown initial conditions, n» 
in the GPM method of Kubicek and Hlavacek (1972). In this example, 
n = (Xh(to).X^(tQ),Ae(to))* the initial values of the three adjoint variables. 
The e = 0 solution would be a starting point. (Note that this means L* = 0.) 

Other parameters are possible; for example, Breakwell (1977) uses 
e' = 1/(2(L/D)„^) for the same problem. This is a coefficient of the induced 
drag term, so that the Hamiltonian would be linear in L for e' =0, producing 
a singular arc. Both formulations will be considered. 

The boundary conditions are 


h(to) = h(j. 

y(to) = Yo(e) 

e{t(,) = eg. 

h(t^) = t^ 

e(t^) = e^. 

\(tf) = 0 


The vector n defines the unspecified initial conditions. 

This problem, then, is in a form to be solved by GPM, starting from £q = 0 
(or £q). The algorithm involves integrating the parametric system 




(n»e) 


iF 

9e 


where 

F^.(x(tf,n>e)) = 0 


is the 


i^ terminal 


boundary condition resulting from an integration of the 


37 


dynamic system to t^ from a given n» and 


r^(nje) = the Jacobian {3F^./9nj} 

In this formulation the integration occurs for all e in where Cj, 

U T T 

is the desired value for e. Note that, at each step in the integration, it is 
necessary to integrate the dynamic system to t^,, in general. (This may be re- 
laxed to every k- step of e if convergence is good.) The differential equations 
arise from applying the implicit function theorem to the system F^=0. 

Ultimately, it is desired to designate e as the singular perturbation para- 
meter. The application of schemes similar to GPM to this problem is not as 
straightforward, because in this role, e = 0 reduces the order of the system. 
However, it is possible to transform the problem so that the scaled time para- 
meter, T“t/e, is the new independent variable in the dynamic system. The 
solution could then proceed from e = 0 as follows: with e = 0, compute the 0 — - 
order solution, including boundary layer solutions at both tg and t^. Boundary 
layer matchings of solutions are necessary to assign values to system parameters 
to insure a stable integration. This 0— -order solution for x(t) could then 
be the first iteration in a continuation process, which would advance e from 0 
to a value small enough so that a first order linearized expansion is adequate. 

We discuss continuation of singular perturbation problems further in 
Appendix 3.1. In addition. Appendix 3.1 discusses the problem of bifurcation 
in continuation. Appendix 3.2 describes a general continuation algorithm which 
we have implemented and which is based on the work of Kubicek (1976) and Keller 
(1977). Appendix 3.3 presents some simple numerical examples of TPBVP's solved 
by continuation. Finally, Appendix 3.4 contains a large bibliography of con- 
tinuation method references. 
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APPENDIX 3.1 

BIFURCATION AND SINGULAR PERTURBATION PHENOMENA IN CONTINUATION 

The basic continuation problem is to solve the equation 

(1.1) 6(x.X) = 0 

for the vector x in terms of the real parameter A. If B is a Banach space and 
G maps BxR into B, then we wish to find the trajectories X-^x(X) in B which 
satisfy 

(1.2) G(x(X),X) = 0 

for all parameter values X. For example, (1.1) might be the nonlinear equation 
for the missing initial values in the TPBVP. In this case, B would be a finite 
dimensional space, B = R^. On the other hand, one might treat the entire trajec- 
tory control and state trajectory as a vector in an infinite dimensional vector 
space B. In that case, (1.1) would be a nonlinear operator equation on the 
Banach space B given by the Euler-Lagrange necessary conditions for the trajec- 
tory optimization problem. It appears that this infinite dimensional point of 
view offers the advantage of numerically more stable solution algorithms. 

The; method of differentiation with respect to a parameter first discussed 
by Davidenko (1953) solves the equation (1.1) for all values of the parameter X 
by differentiating equation (1.2) with respect to X to obtain 

(1.3) H ^ + If (x(X).X) = 0 

Equation (1.1) is first solved at some value of the parameter X, say and 

the solution x(Xq) is used as an initial condition from which to integrate the 

equation (1.3). The integration may proceed as long as the Frechet derivative 

(x(X),X), which is a linear operator from B into B for a given value of x 
ax dx 

and X, is nonsingular and may be inverted to solve for the derivative -rr 

(1.3) . li" (x(X),X) should become singular for some value of X = X^, then a 
bifurcation has occurred at X^, indicating a sudden change in the nature and 
number of solutions of (1.1) in the neighborhood of 
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In Figures 1 and 2 typical bifurcations* are illustrated. In Figure 1 
there 1s one solution of (1.1) for values of A near such that A<A^, and 
there are two solutions of (1.1) for values of A near A^ such that In 

Figure 2 there is no solution of (1.1) for values of A near A^ and such that 
A>A^, but there are two solutions for values of A near A^ such that A<A^. 

The proper treatment of bifurcations is extremely important in applying the 
continuation method to the solution of nonlinear problems for which there may 
be multiple solutions. Such a situation occurs, for example, if the Euler- 
Lagrange necessary conditions allow multiple extremal solutions of the trajectory 
optimization problem. One attractive feature of the continuation method is that 
it will yield all of the multiple solutions--provided that one has an algorithm 
that can handle bifurcations. We have implemented an algorithm due to Keller 
(1977) which handles both bifurcations pictured in Figures 1 and 2. For ex- 
ample, in Figure 1 the algorithm would trace out the left branch until A ap- 


' 9 G 

proached A^ then it would indicate the singularity of (x(A^),Ap and proceed 

to trace the two right branches. In Figure 2 the algorithm would trace out the 

upper left branch until the parameter value reached A.. It would then indicate 

a singularity in ~ (x(A^),A^) and reverse direction of the parameter A to trace 

the lower left branch. This continuation algorithm is described in more detail 

in Appendix 3.2. 

In solving a nonlinear equation G(x) = 0 by the continuation method, one 
first mbeds this problem in a one-parameter family of problems represented by 
(1.1) :'or which the parameter value A=1 gives the original problem and the 
parameter value A = 0 gives a problem with a known solution. If one continues 
al 1 of the solutions of (1.1) at A = 0 to values of A for which A > 0, then one 
is gua 'anteed to find all the solutions of (1.1) at values of A>0, provided 
that one follows all the branches fro m bifur cations . In this way, one finds 
all the solutions of the original nonlinear problem G(x) = 0. In the case of 
an optimization problem, such as the trajectory optimization problem, one ob- 
tains all the extremal trajectories which satisfy the necessary conditions, and 
one may pick out the optimal solution from these extrema. 

In addition to following all branches from bifurcations, it is also impor- 
tant to know all the initial solutions at A = 0 and to continue each of these 


*See note for Figure 2. 
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A 


Figure 1. -Example of Bifurcation Point* 



A 

Figure 2. -Example of Limit Point. 


NOTE: Unhappily, as in the case of the word "singular," "bifurcation" has come 

to have at least two distinct meanings. As opposed to the above definition, sev- 
eral authors, e.g. Keller (1977), designate only the Ai in Fig. 1 as a bifurca- 
tion point, making it a subset of our definition. In this case, A^^ in Fig. 2 
is named a 1 imit point. The cases of Fig. 1 and Fig. 2 are distinguishable 
mathematically, as will be shown in Appendix 3.2. 
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solutions to values of A>0. One usually selects the one-parameter imbedding 
(1.1) so that it yields only one solution x(X) for X in a neighborhood of 0. 
However, it sometimes happens that (1.1) has no solution at A = 0 although it 
has a unique asymptotic behavior as A^O. Such is the situation in problems 
such as our aircraft trajectory optimization problem, in v/hich the parameter A 
represents a singular perturbation near 0. For example, the parameter A is 
denoted by e in the aircraft optimization problem represented by (4.2.46), 
(4.2.47) in Section 4.2. These equations have no solution at A = 0 since it is 
generally impossible to satisfy all the boundary conditions for the reduced or- 
der (A = 0) equations. However, the optimal solution has a well-defined, unique 
asymptotic behavior as A->0. In this case, it is still possible to apply the 
continuation method successfully by starting the process at a A^ > 0 in a neigh- 
borhood of A = 0 and by using the asymptotic approximation of the solution x(Aq) 
as the initial condition in the differential equation (1.3). To do this, one 
may have to make preliminary Newton corrections to obtain a more exact initial 
condition at A^, but for small nonzero values of the parameter A the asymptotic 
approximation is very accurate and only a few Newton corrections are usually 
necessary. The implemented continuation algorithms feature such preliminary 
Newton correction, as well as Newton correction after every prediction step, to 
improve accuracy to a pre-specified tolerance. 

Continuation from a singular perturbation will be effective when the asymp- 
totic approximation is good for A near A = 0 but not very good for A near the 
desired solution at A=l. In particular, this technique will improve the 
asymptotic approximation to the aircraft optimization problem when the total 
range for interception decreases below the lower limit for validity of the 
energy state approximation.* 


*This lower limit is about 160-170 Km for the F-4 example problem. 
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APPENDIX 3.2 

GENERAL CONTINUATION ALGORITHMS 

A series of numerical algorithms for the solution of nonlinear systems of 
the form (A3. 1.1) has been developed for application to such problems as air- 
craft trajectory optimization. Ultimately, it is hoped to combine these al- 
gorithms into one system. In this way, the simplest (and presumably the most 
time/core-efficient) algorithm may be selected initially, switching to more 
complete algorithms as the numerical development of the problem warrants. The 
basis for the developed algorithms lies principally in work by Klopfenstein 
(1961), Kubicek (1973, 1976), Keller (1977) and Rheinboldt (1977). 

Each of the aforementioned researchers realized that effective implementa- 
tion of the parametric system of Davidenko (A3. 1.3) requires adequate treatment 
of the Frechet derivative 

(2.1) F(x(X),^) = II (x(A),A) 

when X is in the neighborhood of some point X^ at which F becomes singular. F 
is singular at either point X^ in Figure 1 of Appendix 3.1 (a "proper" bifur- 
cation point) or Figure 2 (a "limit" point). Rheinboldt and Keller have de- 
veloped a means of extending the continuation through the Figure 1-type bi- 
furcation point. This point differs from the limit point of Figure 2 in a 
mathematical sense, as will be seen below. Limit points are somewhat easier 
to deal with, and so will be discussed first. 

As can be seen in Figure 2 (especially if one imagines x and G to be scalar), 
dx/dX-»-«> as X^X;j^ from Xq. Hence, F must be singular at X^ for the trajectory 
to be meaningful. It is also precisely at x^ where the parameter x loses its 
monotonic property. Because the connection is not coincidental, it seemed 
reasonable to augment the problem by introducing an arclength parameter, say t, 
which is by nature monotonic. The problem is augmented in that the system 
parameter X is now itself a dependent variable in the arclength parameter. Thus, 
if X€ R*^, the augmented system is based on the solution of 


(2.2) G(x(t).X(t)) = 0 
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or 


(2.3) G(y(t)) = 0 
where 

(2.4) y = (x,X), ycR"'^^ 

Since there are now n+1 unknowns, (2.2) must be augmented by the particular 
arclength relationship. Both Klopfenstein and Kubicek use a purely Euclidean 
relationship (called a "normalization") 

(2.5) N^(x,X,t) = 0 = + ... + x^^ + - 1 

where 

( 2 . 6 ) 

The system can now be solved, as follows: 

(2-7) ^-5 = ||x+fx = 0 


Equation (2.7) is the starting point for both methods of Klopfenstein and Kubi 
cek. It may be rewritten 


( 2 . 8 ) 

where 

(2.9) 


Ay = 0 



Equation (2.8) is a system of n equations for the n+1 unknown elements of y. 
Klopfenstein does not distinguish between the x.j and X, considers ^ = 

His method depends on finding an x^^^ for which F remains non-singular. Then, 



Equation (2.5) then determines 

The method of Kubicek is more "robust" in that, while he also considers X 
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and X functionally equivalent in terms of t, the most non-singular nxn sub- 
matrix of A is used at each point for the matrix inversion operation. This 
is achieved by Gaussian elimination with controlled pivoting. By this procedure, 
one of the n+1 columns of A for which A is "most singular" is eliminated, say 
column k. Then X|^-«-»-y|^ plays the role of the parameter A, instead of always 


using Equation (2.7) is then rearranged and solved for the n y^., i^k: 

( 2 . 11 ) 


where 


( 2 . 12 ) 


(2.13) 


~ - F Tv 1 + 

dt 3y^ 


yk = o 


F. = 


y. = -F, 


3Gj 

BGj 

9yi ’ ’ 

• 


3G. 


ayj 

• 




9yi ’ *** ’ 



i k 


3G^ 




k+1 


3G. 


9y 


n+1 


5 

■ ’ 9y 


n+1 


Note that Fj^ is nxn, a square matrix. Equation (2.13) is similar in form 
to (2.10), but more general. 

As before, the y^. are substituted into (2.5) to solve for the final element, 
yj^. The sign ambiguity is handled by selecting a sign at the starting point, 
yQ, t=0, which is consistent with the problem, realizing that t increases 
monotonically. For example, if rudder deflection, fir, is selected as and 
has the value +fir„^^ at the starting point, then the negative sign on y. would 
be selected, so that continuation may proceed into the acceptable range of fir 
val ues. 

The Gaussian elimination procedure works as follows: all of the elements 

of A are scanned for the one with the largest magnitude. This element becomes 
the pivot point for the elementary matrix row operations which are used to zero 


45 



all remaining elements in the column of the pivot element. This column is 
saved and the scanning process begins again, for the remaining columns. This 
process is only done n times, leaving untouched the column whose elements are 
consistently the smallest. The untouched column becomes column k. The remain- 
ing n columns can be rearranged to produce a diagonal matrix, whereupon inversion 
is straightforward. 

Again, the method works well at limit points because, say at in Figure 
2 of Appendix 3.1, X would cease to be the parameter, but would be replaced by 
one of the x^. . It should also be mentioned that provision can be made for in- 
fluencing the choice of by scaling each of the columns of A by a scalar, thus 
reducing their magnitudes. 

Given that the n+1 y^. have somehow been found at a certain point te [0,1], 
the continuation process of Kubicek evolves essentially by integrating the sys- 
tem (2.13) and (2.5) for y(t). The continuation proceeds numerically by a 
predictor-corrector sequence. At the starting point (t=0, X = Xq), the x^ are 
found to the required precision using Newton-Raphson. Then, Adams-Bashforth 
variable order (<4) is used to advance all of the y^. --i.e., x. and X — to a pre- 
dicted value. At this new point, yj^ is found as described above, and Newton 
corrections are made on the y^ , i k, until reasonable convergence is assured. 
Typically, no more that three or four Newton-corrector steps are needed at 
each point. The Newton-corrector formula is 

) 

where y^ * ^ is an n x 1 vector of al 1 y^ , i ^ k. 

The Klopfenstein/Kubicek algorithm was also tested on another project. In 
this application, equilibrium solutions of an aircraft are generated in both 
developed spin and roll departure flight regimes. The continuation system is 
the equations of motion, with the derivatives set to zero. One of the three 
aerosurface controls is chosen as the nominal continuation parameter, and the 
other two are set to fixed values. For spin analysis, an eighth order system 
of equations is needed, because of coupling effects arising from the need to 
employ the full nonlinear expressions. (The roll departure regime requires 
only a fifth order system, but is still highly nonlinear.) However, the con- 
tinuation algorithm described above is able to solve for the equilibrium surfaces 
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(2.14) 


y(j+l) = y(j) 


- F, 




(j) 



quite readily. This algorithm is also the basis for some of the two-point 
boundary value problem (TPBVP) continuation examples which are described in 
Appendix 3.3. For TPBVP applications, the system (2.5), (2.13) operates in 
function space rather than Euclidean n-space. The significant difference is 
that, for every continuation step ("outer" integration), a full solution must 
be generated between the two boundary points by quadrature ("inner" integration) 
While the function-space process is thus considerably more complicated, the 
continuation itself is non-iterative; that is, the continuation parameter moves 
directly from its initial to its final value. There is no retracing or repe- 
tition of this outer loop quadrature. 

For flight trajectory optimization problems, the TPBVP system consists of 
adjoint, or influence, functions whose initial conditions are unknovm in general 
These, then, become the variables x in (2.2) and A, the continuation parameter, 
is usually chosen to be a physical parameter such as (L/D)__^ or air density. 
This choice enables a reasonably simple solution to be found at A = 0. See 
Appendix 3.3 for a discussion of continuation methods applied to the solution 
of TPBVP 's. 

At this point, it is probably easy to appreciate that it is quite often 
convenient to evaluate the elements of A numerically, rather than to pre-compute 
analytic expressions. This is particularly the case when simulating flight 
vehicle trajectories, because such expressions involve terms containing tabular 
data. Analytic expressions for such terms are typically very complicated. 
Therefore, a numerical differentiation algorithm has been developed and is now 
a part of our basic continuation algorithm. It may be invoked at the user’s 
option. 

The algorithm which thus far has worked most efficiently is the following: 
to compute fit a cubic spline to G(x) at the following five 

points (knots): 

(2,15) G(x+iA), i = 0,±l,±2 

Once this fit is made, the slope of the spline function at Xq becomes the ap- 
proximation to F(Xq) = Fq. In genera], G may be a vector and x a scalar. If 
xeR^, the extension is obvious; one merely performs the operation n times for 
each X.. 

\J 
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It is important to use a value for A which is small enough to adequately 
represent yet not so small that numerical difficulties result. Thus, we 
require* 

(2.16) |G(x) - G(x + A)| 0(e) 


v/here e = 0.0001 has been shown to produce good results. 

The algorithm for selecting A is as follov/s: 

1) initial ize i = 1 

2) set A^ = G (this is exact if F(xq) = 1) 

3) compute G. = |G(Xp) - G(Xq + A^ ) | 

4) lOe, go to 7 

5) if not, set 

Note that 6^. = 0 where the slope is infinite 

6) Set i = i+1 and go to 3 

7) set A = A. and do the spline fit 

If a solution for A is not obtained after five iterations, a v^arning is printed 
and A is set to A^. 

As mentioned above, the Kubicek algorithm is unable to solve automatically 
for all of the branches which emanate from the bifurcation point shown in 

Figure 1 of Appendix 3.1, We have implemented, therefore, a method based on the 
v/ork of Rheinboldt (1977) and Keller (1977) which can not only continue the 
original branch accurately past but also accurately evaluate X^ and the 
slope of the "secondary" branch at X^. With this slope, continuation along the 
secondary branch can proceed as usual . 

Keller's algorithm begins with the system (2,2) augmented by the "pseu Jo- 
arc length" normalization 


(Z.18) N3 = exQ^(x-Xp) + (l-e)A(i-AQ) - (t-tp) = 0 


where t is the arclength parameter, 8 is a constant selected such that 0<6<1, 
and (XqjAq) = yg are the values of the unknowns and the parameter at tg. The 

*If G f r", (2.16) is modified: 

(2.17) i iG.(x) - G (x + A)| = 0(e) 

" i=l ^ ^ 
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augmented system 1s novr. 


(2.19) 


H(y,t) 


6 


'G(y) 

Nq{y,t) 


= 0 


where 


(2.20) y = [xl ’ y^ 

Keller has shov/n that, using as given by (2.18), the quantity 


(2.21) B ^ ly 


is nonsingular if and only if 

a) nonsingular; or 
3G 

b) lx where R(0 denotes range space. 

Case (b) corresponds to a limit point. At such a point (e.g.. Figure 2 of Ap- 
pendix 3.1) there is no intersection of branches, but dx/dA-^<». However, solu- 
tion of the augmented system (2.19) continues normally. 

If neither condition (a) or (b) holds at some point then B is singular 
and Xj is the type of bifurcation point shown in Figure 1, v^here two or more 
branches intersect. As mentioned above, is skipped over, in order to continue 
along the initial branch. Any predictor-corrector method will suffice for this. 
Continuing the solution along the second branch then proceeds as follows*: 

1) compute 3G/3X (singular) at the bifurcation point, and get an approxi- 
mation for dx/dt and dX/dt. 

2) compute the eigenvalues and eigenvectors of 9G/9x. Call the eigenvector 

corresponding to the zero eigenvalue 4)g, and the others • 

3) compute (|)j^ in the following: 

(2.22) (3G/9x)(})j^ + 36/3X = 0 

where <|)j^ is a linear combination of the vectors ^ 

vectors (|)q and define the plane in v/hich a solution for 9x/9t can 

appear. 

4) solve for the coefficients Tq and in: 

*This is Keller's (1977) "Method II." 
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(2.23) 


^0 = 

" ■'^0 

6) novi we seek a solution of the form 

(2.25) = Xg + e[’fo'!'o + ‘fN<t>N] ■*■ 

" ^0 e^N 

where Xq,Xq is the solution at the bifurcation point. 

7) final ly» starting with V=0, n=0, use Newton's method to find a solu- 
tion to the equations: 

(2.26a) &(x\x^) = 0 

(2.26b) + t^4,*)V + = 0 

The c terms in (2.25) move the possible solution av/ay from the first branch, and 
equation (2.26b) insures that the free variables V and n do not bring the solu- 
tion back. Thus, any solution found must be on the second branch of the solution 
curve, and can be used as a starting point for a normal continuation solution. 

This method has been implemented and run on test problems. Because of its 
complexity, especially in the use of TP8VP continuation algorithms, much effort 
required to select numerical parameters and quadrature techniques which will 
efficiently deal with the system at hand. 
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APPENDIX 3,3 

APPLICATION OF CONTINUATION ALGORITHM TO TPBVP: NUMERICAL EXAMPLES 

The algorithm described in Appendix 3.2 was applied to solving three ex- 
ample TPBVP ‘s by the shooting method. The theory behind the shooting method is 
presented in Subsection 3.3.3 and we refer the reader there for details of the 
basic method. Note that the present method of solving the TPBVP by shooting 
uses both the Davidenko method of differentiation with respect to a parameter 
described by Kubicek and Hlavacek (1972a, 1972b, 1973) and Kubicek (1976) and 
the Newton-Raphson method described by Roberts and Shipman (1967, 1968) to- 
gether in one algorithm. 

Example #1 

This is a nonlinear second order two- point boundary value problem which we 
solved by continuation with respect to the final time as described in Roberts 
and Shipman (1967, 1968, 1972). The problem was to solve 

dUj 

dt~ ° “2 

du^ 

(3.1) ^ In Uj 

Uj(0) = 1, = e 

for several given values of the final time T, The problem was solved numerically 
by continuation with respect to the final time T, starting at T = 0 where the 
problem reduces to an initial value problem with u^(0) = 1 and U 2 ( 0 ) = e. The 
problem in (3.1) was chosen to be the nonlinear transformation of the linear 
probl em 

dx, 

dt~ " ^2 
dx^ 
dt” "" 

Xi(0) = 0, X2(T) = 1, 
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via the transformation and so that we would have an exact solu- 

tion to the problem (3.1). If b(T) denotes the missing initial condition 
U2{0)=b(T) in (3.1) for the problem with final time T given, then the function 
b is given by the formula 

(3.3) b(T) = 8xp(2/(e^ + e'^)) 

The numerical continuation solution agreed with (3.3) exactly in all displayed 
decimal digits. 

Example #2 

This is also a nonlinear, second order TPBVP but the final time is fixed 
and an internal parameter is used for continuation. This example represents the 
equilibrium equation for the heat distribution in a rod of length x=l. The 
problem is solved by continuation in Wasserstrom (1973) which is an excellent 
review of continuation methods in general. The equations for the problem are 


yi(0) = y^d) = 0 


where X is the continuation parameter used. We applied continuation to the 
shooting method to calculate numerically the missing initial condition 
y2(0) = b(x) for values of the parameter X between 0 and 1. Note that for X = 0 
the problem has a trivial solution, y(x) = 0 for all x. 

The exact solution of (3.4) is given by the formula 

(3.5) yM = ln([in^m]sec^(|[x-y)) 

where m must satisfy the implicit equation 
2 

I f \ m 2m , 

(3.6) 2X 4 ^ 

The missing initial condition b(X) is given in terms of m and X as 
52 



(3.7) b(X) = -m tan(n/4) = -/2A-m 


Again, the numerical continuation solution agreed exactly with (3.7) in all 
displayed decimal digits. 

Note that our continuation algorithm used numerical differentiation to 
calculate the partial derivatives necessary for the shooting method instead 
of integrating four extra variational differential equations as Wasserstrom 
does in his paper. The numerical differentiation simplifies the prograirsning 
requirements and v/ill reduce the computation required in large order problems. 



The last example is a trajectory optimization problem for a simplified 
missile interceptor taken from the paper of Schneider and Reddy (1974). Al- 
though the aerodynamic model is extremely simplified, this problem is a good 
example with which to test out the continuation algorithm before attempting the 
aircraft equations. The state equations in this problem are nonlinear and 
fourth order, resulting in a nonlinear, eighth order TPBVP. The nonlinear, 
fourth order state equations are 


dt ' ^3 
dx2 

dt“ " ^4 

(3.8) 

^ - -«3ve-Vhs . 
dx 

^ = -ex^ve"^2/hs + _ 9^31 

where v = V(x^^ + x^^) . 

In the equations (3.8), represents a horizontal range variable and 
is the height variable. The variables and are the velocities correspond- 
ing to x-y and X 2 respectively. The variables Uj and are the controls 
(thrusts) for the problem. The parameter h is a height scale which is taken 
as 6705.6 m in this problem- The parameter e is equal to ^ where is the 
sea level air density and b is the Inverse ballistic coefficient for the missile 
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interceptor modeled by (3.8). The optimality criterion for this problem is 
the simplest quadratic criterion > 


(3.9) 


f (u,^ + U2^)dt 


where T is the fixed final time. The Euler-Lagrange necessary conditions give 
the following nonlinear, eighth order TPBVP: 


dt ^3 

dx^ 

7T " 


(3.10) 


dx. 

dt 

dx^ 

dT 

dX( 

dt 


'GX^ve 




- 9.81 


= 0 


dt 


ev 

h. 


e-X2/hs^ 


V4 


dt 

dt 


"^5 

■^6 




where Xj(0), X2(0), x^lO), x^(0), Xj(T), X 2 (T), are given and X 7 (T) = Xg(T) = 0. 

As in (3.8), v = V(^ 3 ^ + . In equation (3.10), the variables Xg, Xg, Xy and 

Xg are the adjoint variables corresponding to x^, x^, Xg and respecti vely. 

We solved (3.10) by shooting for the missing adjoint initial conditions Xg(0), 
Xg(0), Xy(0) and Xg(0). The parameter e (proportional to the inverse ballistic 
coefficient) was used for continuation from e = 0. Note that the g = 0 solution 
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corresponds to an interception problem with no air density — i.e., a vacuum 

solution. Figures 1, 2, 3 and 4 show the range-height trajectories for four 

4 -2 

values of the ballistic coefficient ranging from to 3.1266x 10 Nm . This 

corresponds to continuing the parameter e from 0 {the vacuum solution) to 
-4 -1 

1.9220x 10 m . The final time T is 20 sec and the given initial and final 
conditions are X2(0) = X3(0) = x^(0) = 0, x^(0) = -1.524 x lo^ m, x^(T) = 0, 

X 2 (T) = 4.409 X 10 m. These initial and final conditions correspond to the con- 
ditions for cases 5 and 6 in Schneider and Reddy (1974). Schneider and Reddy 
presented the trajectory of case 6 which corresponds to our Figure 4, although 
the ballistic coefficient for case 6 is about two times that of Figure 4. 

Figure 3 is very close to case 5, although the ballistic coefficient for Fig- 

4 -2 

ure 3 is 4.6923 xio Nm which is slightly less than the coefficient for case 5, 

4-2 R 

namely 5.788xio Nm” . The cost associated with Figure 3 (A= 3.542 xio ) is 

Q 

correspondingly slightly larger than the cost in case 5 (A=3.415xio ). 

As the parameter e and the final time T increase, the equations (3.10) 
become more sensitive to slight changes in the initial values of the adjoint 
variables. This is partly due to the increased nonlinearity of the problem, 
but the main trouble comes from the forward integration of the adjoint equations. 
The adjoint equations are unstable in the forward direction and the differential 
equations (3.10) will become infinite in a finite amount of time. One solution 
to the problem which still maintains the shooting method is to use double pre- 
cision instead of single precision accuracy in computations. A1 ternati vely, 
one can use one of the function space methods such as quasilinearization as 
presented in Roberts and Shipman (1968a, 1972) or the back-and- forth shooting 
method of Orava and Lautala (1976, 1977). 
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SCHNEIDER-REDDY EXAMPLE 

Figure 1 
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SCHNEIDER-REDDY EXAMPLE 

Figure 2 
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INTERCEPTOR BALLISTIC COEFFICIENT (1/b) = 4.6923 x 10^ Nra' 
A = 3.542x 10® 



Y(1 ) 

RANGE (units OF 6705.6 M) 
SCHNEIDER-REDDY EXAMPLE 


Figure 3 


HEIGHT (units OF 6705.6 m) 
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CHAPTER 4 

SPT APPROXIMATION OF OPTIMAL CONTROL LAWS 


4.1 Introduction 

For the purposes of on-board control of an aircraft we are not interested 
so much in obtaining optimal control trajectories as in obtaining optimal control 
laws. The essential difference between these two ways of specifying cdntrol is 
that the optimal control trajectory specifies the control to use as a function 
of time while the optimal control law specifies the control to use as a function 
of state. In principle, the control law is more valuable because it can provide 
feedback correction to perturbations in the optimal trajectory caused by uncer- 
tain environmental factors, e.g., wind gusts, or by errors in the mathematical 
model of the system, e.g., modeling the aircraft as a point mass rather than 
as a rigid body. Computationally, however, control trajectories are much easier 
to compute than optimal control laws because one need only calculate a function 
of one time variable rather than a function of several state variables. Indeed, 
even if we could calculate the optimal control law exactly off-line, the storage 
requirements would prohibit us from using this as an on-board control for any 
but the most modest sized problems. Suppose n is the number of states for the 
system, m is the; humber of controls, and N is the number of discrete values we 
divide any one of the n state variables into for storage. Then to store the 
control law requires storing an n-dimensional array with n x elements. In our 
problem we have n=6 and m=3. Even if we were willing to accept the crude 
approximation of the control law resulting from using only N = 10 discrete values 
for each of the n = 6 state variables, we would still have to store a 6-dimensional 
array with 30xio elements. Thus, the well-known "curse of dimensionality" 
forces us to use approximations for both the computation and the implementation 
of optimal control laws. 

In this chapter we present some general methods for approximating optimal 
control laws by means of singular perturbation theory (SPT) techniques. In 
the next chapter we apply these general methods to the minimum time control of 
the six-state aircraft model described in Chapter 2. A brief summary of the 
contents of the present chapter follows. 
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Section 4.2 describes the basic calculation of the SPT approximation to 
the optimal control law for a general dynamic system. To clarify this cal- 
culation we discuss in detail the case of two time scales. The general case of 
several time scales is discussed more briefly — we present the results and leave 
the details of the derivation to an appendix. Section 4.3 discusses the compu- 
tational difficulties of the SPT calculation described in Section 4.2 and it 
also discusses alternative strategies for overcoming these difficulties in order 
to obtain efficient real-time algorithms. In particular, we discuss in detail 
the use of individual time scales for each individual state variable, lineari- 
zation of the SPT calculation outside of regions of rapid variation and the use 
of suboptimal solutions for the slow time scales. Finally, Section 4.4 discusses 
the problem of the validity of the SPT approximation in various regions of state 
space. For example, by properly scaling the aircraft equations we obtain a 
natural singular perturbation parameter that is inversely proportional to total 
range in the interception problem. In this case, the SPT approximation is better 
the farther the state is from the terminal target states. Near the terminal 
target the SPT approximation breaks down dramatically and implementation of the 
SPT control law leads to a destabilizing feedback law. Thus, in the vicinity 
of the target state one must switch from SPT to some other method of approxima- 
ting the control law. 


4.2 Calculation of the SPT Approximation 

4.2.1 Control Trajectories, Control Laws and the SPT Approximation 

In Section 4.2 we will discuss the SPT approximation of optimal control 
laws of the following generic, autonomous control problem. Let x represent an 
n-dimensional state vector trajectory and let u represent an m-dimensional 
control vector trajectory. The dynamic equations for the system are: 

(2.1) ^=f(x,u) 

Suppose that u is a given control input trajectory and let x be the corresponding 
state trajectory obtained from solving the differential equation (2.1) over the 
time interval [tg,«>) with the initial condition 
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(2.2) x(tQ> = C 


The cost functional J(u) corresponding to the given input u is defined as the 
integral cost criterion 


(2.3) 



L(x,u) dt 


where t^ represents the final time at which we require the state to lie in the 
terminal target set. For simplicity we assume that the terminal target set con 
sists of the single state oj so that the terminal condition for the problem is 


(2.4) x(t.p) = w 


In this generic control problem we will allow the terminal time t^ to be either 
free or fixed and we will specify which is the case unless the results are the 
same for either case. 

The optimal control problem is to find a control input u* such that (2.3) 
is minimized by u = u* subject to the condition that (2.1), (2.2) and (2.4) also 
hold. For a given fixed terminal time and a given fixed terminal condition 

(2.4) , the optimal input u* is a function of the initial condition the ini- 
tial time tg and the current time t where tQ^t<t^. That is, we have 

(2.5) u* = u(t,tQ,?) 

By fixing the initial time tg and the initial state‘s, we can consider the op- 
timal control in (2.5) as a function of the current time t only, i.e., as a 
control trajectory . On the other hand, if we let the initial time tQ and the 
initial state 5 vary and if we set t-tg, then we can consider the optimal 
control in (2.5) as a function of time and state only, i.e., as a control law . 
Here we are interested in obtaining an SPT approximation of the control law 
and we will do this as follows. Using the SPT methods of trajectory optimiza- 
tion, we will approximate the optimal control function in (2.5) for arbitrary 
initial times and initial states. We will then set t = tg in (2.5) to obtain 
the optimal control u* to use at time tg with the corresponding state Since 
tg and C were arbitrary, this procedure will give us an approximate control 
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law to use at any time for any state for which the SPT approximation is valid. 
In Section 4.4 we will see that the SPT approximation is not valid for certain 
values of 

Suppose that the control problem in (2.1) -(2.4) has only two time scales, 
fast and slow, and suppose that x= (x ,x ) where x represents the fast com- 
ponents of x and x represents the slow components of x. Then the singular 
perturbation approximation to (2.1) is 

(2.6) ^=f^(x,u) 

dx^ 2 

(2-7) 6^ = r(x,u) 

where e is p parameter varying from 0 to 1 and f^ is the component of f corres- 
ponding to the x^ component of x, i = 1,2. For each e the optimal control 
problem with (2. 6), (2. 7), (2. 3), (2. 4) results in a TPBVP similar to the unper- 
turbed (e=l) problem. In terms of the Hamiltonian H defined 

(2.8) H(x,X,u) = L(x,u) + X^f(x,u) 

1 2 

where X=(X ,X ), we have the perturbed equations 


(2.9) 

dt 

9H 

1 

3x^ 

(2.10) 

dX^ _ 
dt 

m 

dx' 


which correspond to (2.6) and (2.7). The SPT approximation of the control tra- 
jectory results from approximating the TPBVP given by (2.6), (2.7), (2.9) and 
(2.10). Typically, this approximation has three parts: (1) an initial boundary 

layer at the initial time tQ, (2) a reduced order, or outer, solution for t such 
that tQ<t<t^, and (3) a terminal boundary layer at the terminal time t^. 

These different parts are determined as follows. The initial boundary layer is 
given by solving 
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(2.11) 



(2.12) 



(2.13) 

- 0 


(2.14) 

^ 3x 

t-t 

L Lq 

where t is the fast time scale given by t ~ ^ . The 

this TPBVP are the same as the initial conditions (2.2) 
(e=l) problem. That is. 

initial conditions for 
for the unperturbed 

(2.15) 

S^(0) = 5^ 


(2.16) 

5^(0) = 5^ 


The terminal conditions are not those for the original 
we have the asymptotic terminal conditions 

problem, however. Instead 

(2.17) 

1 im x^(t) = x^(0) 


(2.18) 

lim x^{t) = x^(0) 


(2.19) 

11s = x^(o) 


(2.20) 



where x 

x^, X^, X^ are the states and ad joints for the reduced order solution. 

That 1s 

, x^, x^, X^, X^ satisfy the TPBVP 


(2.21) 

^=fl(x.0) 



( 2 . 22 ) 


0 = f^(x.u) 
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_2 

(2.23) ^ = - -^.(x.X.u) 

(2.24) 0 = - -^-(x.A.u) 

3x 

The initial condition is 

\ 

(2.25) x^(0) = : 

and the final condition is given by 

(2.26) xHtf) = u)^ 

-2 

Note that the reduced order solution treats x as a pseudocontrol which satisfies 
the algebraic equation (2.22) rather than the original differential equation. 
Finally, the terminal boundary layer is given by solving 

dx^ 

(2.27) ^=0 

(2.28) gl=f2(x,a) 

(2.29) ^.0 
da 

~2 

(2.30) 

t-to 

where a is the fast time scale o = . The initial conditions for this TPBVP 

c 

are the same as the final conditions (2.4) of the original problem. That is, 

(2.31) x^(0) = J 

(2.32) x^(0) = J- 

The terminal conditions are given asymptotically in terms of the reduced order 
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solution as 


(2.33) 

(2.34) 

(2.35) A^(a) = A^(t^) 

(2.36) A^(a) = A^(t^) 

After solving (2. 11) - (2.36) we obtain three control trajectories: the 

initial boundary layer trajectory u(t), the reduced order trajectory u(t) and 
the terminal boundary layer u(a). The SPT approximation which is uniformly 
valid over the whole interval tQ<t<t^ is given by 

(2.37) u*(t,tQ,C) - [u(t) - u(0)] + [u(a) - u(t^)] 


JiOL x^(®) = x^(tf) 


t-tQ t-t^ 

where 't=— ^ — and o = ~ — . Note that we have suppressed the dependence of u, 
u and u on t^ and ^ for clarity. Here we have only sketched the SPT method of 
approximating solution trajectories of the TPBVP one obtains for an optimal 
control problem. In Subsection 4.2.2 we will present this approximation in 
more detail for the problem of obtaining control law approximations. The 
reader should refer to Kelley (1973) for more details about obtaining optimal 
trajectory approximations by means of SPT. 

Because we only wish to obtain the approximation of the time function 
t-^ u*(t,tg,^) at the initial time t = tg, we only need to obtain the initial 
boundary layer terms in the SPT approximation. However, to obtain this term we 
need to calculate the reduced order approximation first. Nevertheless, we do 
not need to aonstder the terminal boundary layer at all to obtain the lowest 
order approxi?nation. To obtain the higher order terms of the initial boundary 
layer we would have to calculate higher order terms of both the reduced order 
approximation and the terminal boundary layer approximation. The following 
discussion will help clarify this situation. 

Equation (2.37) represents the first term in an asymptotic approximation 
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for thef optimal control trajectory u* as a power series in e. The full asymp- 
totic series has the form 

(2.38) u* - \(tQ)] + Oj^{t) + [u,^(a) - u,^(t^)]}e*^ 

h 

Where U|, denotes the k— order term of the initial boundary layer approximation, 

^ ■ f"h 

Uj^ denotes the k— order term of the reduced order solution and denotes the k — 
order term of the terminal boundary layer approximation. Note that we have 
suppressed the ^ dependence of the functions Uj^, Uj^, Uj^ in order to maintain 
notatfonal clarity. 

The boundary layer terms Uj^ and Uj^ in (2.38) must have the following 
asymptotic properti es : 

(2*39) 14m Q|^(t) = U|^(tQ) 

(2.40) = U|^(tf) 


Hot'eovei', for reasonably well-behaved problems the convergence in (2.39) and 

(2.40) will be exponential. That is, for some positive constants we will 

havfe 


(2.41) |uj^(t) - Uj^(tp)l < e 

(2i42) lii|^(a) - Uj^(t^) | < e 


for sufficiently large t and |a|. Using the information (2.42) and evaluating 
(Z.38) at t=tg. we obtain 


(2.43) 


u* 


T, {uj^(O) 
k=Q 


k 

c 


+ Ove 




-bk(Vto)/c 



where, as usual, the order notation 0(c})(e)) denotes a function which decreases 
as fast as (p(e) as e->0. The SPT approximation in (2.38) is only an asymptotic 
approximation as e^O. Moreover, in (2.43) the exponential terms 
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^/-bk(tf-to)/0 

the exponential 


)e^ decrease exponentially faster than the terms Uj^(0)e*^. Thus, 
terms are negligible in the asymptotic expansion, and we have 


(2.44) u* = E QjO)e‘' 
k=0 


From (2.44) we see that in order to obtain the k— order approximation for.u* 

•f" ki 

we need only find the k— order initial boundary layer approximation 


k . 

(2.45) E Ui(0)e'^ 
j=0 J 

However, note that to obtain the term Uq in (2.45) we need to calculate the 

reduced order term first. To obtain the next term u. we must calculate u^ 

- ^ th ^ 

and u^ as well. In general, to obtain the k— order initial boundary layer 

approximation we have to calculate the reduced order approximation up to the 

th *f"h 

k— order and the terminal boundary layer approximation up to the (k-1)— order. 

“h h 

The important point to note is that to find the 0— order SPT approximation, the 
one which we will be using to approximate the aircraft control law in the succeed 
ing chapters, requires only the calculation of the initial boundary layer and 
the reduced order approximation of the 0— order — the terminal boundary layer 
approximation does not enter the calculation at all. 

The expression (2.43) indicates when the SPT approximation breaks down; 
namely, when the initial time t^ and the final time t^ are very close together 
so that (tQ-t^)/e is small and the exponential terminal boundary layer terms in 
(2.43) cannot be neglected. In this case, the SPT approximation of u* in terms 
of initial boundary layer, reduced order and terminal boundary layer approxi- 
mations is invalid and the asymptotic approximation in (2.38) is incorrect. In 
particular, one cannot improve (2.44) by including the terminal boundary layer 
terms in (2.38). We discuss this situation further and illustrate it with 
simple examples in Section 4.4. In the next subsection we discuss the details 
of the calculation of the k = 0 term of (2.44) for the two time scale case. 

4.2.2 Calculation of th e S PT Approximatio n: T wo Tim e Case 

Although the aircraft example considered in this project is treated as"a 
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system with more than two time scales, the notational difficulties of presenting 
the multi-time scale case obscure the basic simplicity of the SPT calculation 
of the approximate optimal control. Therefore, to illustrate the basic method 
we will discuss the case of a two time scale dynamical system in some detail in 
this subsection. In the next subsection, 4.2.3, we will present and discuss the 
results for the multi -time case, but we leave the details of the calculation to 
Appendix 4.1 at the end of Chapter 4. 

Consider the dynamic system 

(2.46) ^ = f(x,y,u) 

(2.47) = g(x.y.u) 

where xe X, yeY and ue U. We assume that X, Y, U are vector spaces which are 
not necessarily one-dimensional. The cost criterion for the control problem is 

(2.48) J(u) - r\(x,y,u) dt 

where t^ is a terminal time which may be either fixed or free. The terminal 
conditions for the problem are 

(2.49) x(t^) = x^ 

(2.50) y(tf) = y^ 

and the initial conditions for the problem are 

(2.51) x(0) = C 

(2.52) y(0) = n 

To find the first term of the asymptotic series (2.44) we need to obtain 
the initial boundary layer approximation, but as we indicated in the previous 
section, we must first calculate the reduced order approximation. The reduced 
order approximation to the system (2.46) - (2,52) is obtained by setting the 
perturbation parameter e = 0 in (2.47) and omitting the initial condition (2.50) 
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and the terminal condition (2.52). If the quantities x, y> u denote the reduced 
order approximation, then we obtain these quantities by solving the problem with 
the dynamic system 

(2.53) || = f(x.y.u) 

where y and u are both considered as controls for the state x. In addition, we 
require that the following equality constraint be satisfied: 

(2.54) 0 = g(x,y,u) 

The cost criterion for the reduced order problem is the same as for the original 
problem (2.48), namely 

(2.55) 3(y,u) ^ L(x,y,u) dt 

"^0 

In the reduced order problem only the initial and terminal conditions for x are 
retained. Thus, we have the initial condition 

(2.56) x(0) = ^ 

and the terminal condition 

(2.57) x(t^) = x^ 

To solve the reduced order problem we define the reduced order Hamiltonian 

(2.58) H(x,y,u.X^) = L(x,y,u) + X^f(x,y,u) 

in which X denotes a reduced order adjoint variable for the state x. The 
minimum principle implies that the optimal controls y* and u* are chosen so that 
on the optimal trajectory we have 

(2.59) H(x(t) ,y*(t),u*(t) ,A (t)) = min H(x(t),y,u,X (t)) 

^ y,u 

Provided that we can solve the minimization problem in (2.27) for the controls 
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y* and u* in terms of x and X , we have the following two-point boundary value 
problem (TPBVP) for the x and A trajectories: 

(2.60) f(x.y*(x,y,u*(x,y) 

( 2 . 61 ) ^ - H 

The boundary conditions for this TPBVP are the initial and terminal conditions 
on X given in (2.56) and (2.57). If the terminal time t^ is free, then there 
is an additional condition on the Hamiltonian, namely, 

(2.62) H(x(t),y(t),u(t),Xj^(t)) = 0 
for all times t such that 0<t<t^. 

In order to obtain the initial boundary layer we need to have the optimal 
reduced order adjoint X* and the optimal reduced order pseudocontrol y* at the 
initial time t = 0. Assuming that the terminal state (2.57) is fixed and that 
the initial condition (2.56) is variable, then the adjoint X* and the pseudo- 

J\ 

control y* at the initial time t=0 depend on the initial state ^ and the final 
time t^ as follows: 

(2.63) X* = X*(5.tf) 

(2.64) y* = y*(S.tf) 

If the final time is free, then the adjoint and control at the initial time t = 0 
depend only on the initial state C- Having obtained (2.63) and (2.64), we can 
now formulate the initial boundary layer problem. 

The initial boundary layer problem is obtained from (2.46) and (2.47) by 
transforming these equations to the fast time scale T = t/e. Thus, we have 

(2.65) ^=ef(x,y,u) 

(2.66) ^=g(X,5;,(j) 
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where y» u denote the initial boundary layer approximation to the quantities 
X, y and u. Note that in this approximation both x and y are states and u is 
the only control. Since we are only interested in the 0— order approximation, 
we set e = 0 in (2.65) to obtain 

(2.67) ^=0 

(2.68) ^=g(S,J,0) 

Thus, the state x is constant given by the initial condition (2.51), 

(2.69) x(t) = ^ 

for all T>0. Likewise, the adjoint corresponding to x is constant and is 
given by 

(2.70) X^(t) = \(?»V 

for all T>0. Using (2.69) and (2.70), we may rewrite the initial boundary 
layer problem as an infinite horizon control problem for y and u alone. The 
dynamic system for this control problem is 

(2.71) ^=g(^,y,u) 

with the initial condition 

(2.72) y(0) = n 

which is derived from (2.52). The terminal condition for the problem is given 
as an asymptotic limit, namely 

(2.73) :[im P(t) = y(5.tf) 

This asymptotic terminal condition is the requirement that the initial boundary 
layer match up with the reduced order approximation. The infinite horizon cost 
criterion for the initial boundary layer problem is given in terms of the cost 
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function L and the reduced order adjoint as 

(2.74) J(ti) = 

The optimal control trajectory u* for the initial boundary layer problem depends 
directly on the initial condition (2.72) and the fast time t, and indirectly on 
the initial condition ^ and the final time t^ for the reduced order problem. 
Thus, 

(2.75) u* = u*(T,n,?,tf) 

As we saw in (2.44), we only need (2.75) at t = 0 for the SPT approximate optimal 
control law. Thus, the approximate optimal feedback control law u* is 

(2.76) u* - u*(0,n,?,t^) 

The approximation (2.76) expresses the control u* as a feedback function of the 
current state (C>il) and the time-to-go t^. 

In terms of the preceding discussion the calculation of the SPT approximate 
optimal control law can be summarized as follows (also see Figure 4.2.1): for 

a given state (^,n) and time-to-to t^ 

(i) calculate the optimal control adjoint and the pseudocontrol law for the 
reduced order problem in (2.53) through (2.57). Note that the adjoint 
and the pseudocontrol for the reduced order problem will depend on the 
current x state ^ and the time-to-go t^. 

(ii) Using the adjoint and pseudocontrol from the reduced order approximation 
to define the cost and the asymptotic terminal condition for the initial 
boundary layer problem, calculate the optimal control law for the initial 
boundary layer problem (2.71) through (2.74). Note that this solution 
will depend on the current y state n directly, and indirectly on the 
current x state ^ and time-to-go t^. The current x state and the time-to- 
go enter the initial boundary layer problem only in the asymptotic terminal 
condition (2.73) and the cost functional (2.74). 

Before turning to the general multi -time scale case, let us note that without 
further assumptions the SPT approximation is no easier to solve than the 


J [L(^,y,u) + A^(C*t^)f(S,y,u)] dT 
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adjoint and 

pseudocontrol y*(^,t 4 r) 


Figure 4.2.1 .-General SPT Algorithm: Two Time Scales 
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original control problem in (2.46) through (2.52). Although the reduced order 
problem in x is a smaller dimensional problem than the original problem in x and 
y, the initial boundary layer problem in y could be as difficult to solve as 
the original problem due to the general cost functional (2.74) and the general 
terminal condition (2.73) which depend on the current x state We will discuss 
this point in more detail in Section 4.3. 

4.2.3 Calculation of the SPT Approximation: Mu l t j:- Ji m^ Case 

The calculation of the SPT approximation to the optimal control law in the 
multi-time case is conceptually the same as in the two-time case. Therefore, 
we will present only the basic algorithm for the calculation in this section 
and leave the details of the formulation and calculation to Appendix 4.1 at the 
end of Chapter 4. Consider the generic, autonomous control system introduced 
in Subsection 4.2.1 in equations (2.1) through (2.4). Let us suppose that the 
state vector x has been decomposed into lower dimensional components x^ of 
dimension n. for i = 0,l,...,r. If x has dimension n, then 0<r<n-l and for 
each i we must have l<n^<n. Each component corresponds to a separate time 
scale with Xq representing the slowest time scale and x^, representing the fastest 
time scale. The time scales are arranged in order so that for each i the time 
scale for x^-^^ is faster than the time scale for x^. . We now present the al- 
gorithm for calculating the SPT approximation to the optimal control in this 
case of r time scales. The reader may refer to Appendix 4.1 for details of the 
singular perturbation theory formulation and derivation. 

To facilitate our presentation we introduce the following notation. A 
superscript i will always refer to a vector component corresponding to the x^ 
components of the x vector. Thus, we write f^ for the components of the vector 
function f in (2.1). On the other hand, the subscript i will refer to quantities 

■Mn 

corresponding to the i— boundary layer problem. Thus, we define the pseudo- 
control u^. as 

(2.77) u.j = x'",u) 

for i = 0,1,. . . ,r-l and 

(2.78) = u. 
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Similarly, we define corresponding to the first i+1 initial conditions, 

(2.79) 5,- = (5°.....?’) ■ 

for i = 0,1,. . . ,r. 

The vector functions g. of (x,u) are defined so that 

(2.80) g.j(x,u) = (f^'^^(x,u),. . . ,f’^(x,u)) 

for i = 0,1 r-1 and so that 

(2.81) g^(x,u) =0 

The function g^+;^» which is trivially 0, is used only for notational cpnvenience, 

A 1 gorithm for SPT Approximate Optimal Control Law* 

The control to use with current state ? and time-to-go t^ is calculated as 
follows: 

(i) Solve the reduced order problem: 

(2.82) ^ = f°(x°.Uo) 

(2.83) 0 = go(x°,U(,) 

(2.84) x°(0) = 5° 

(2.85) x°(tf) = lP 

where w is the terminal condition (2.4) for the original problem, 

(2.86) Jq(Uq) = L(x°,Ug) dt 

■'o 

where L is the cost criterion (2.3) in the original problem. 

(ii) Let x^ denote the first component of the pseudocontrol Uq, and let 

*Also see Figure 4.2.2. 
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Xq, 


Solve initial boundary layer control problem 
#=fnCi.i,x\up 
0 = 

x^(0) = x^(”) = x^ 

•-1 = L,--l " 

'Jo h(5,_i.x’.u^) dt 


i* ^ 


Approximate SPT optimal control 
u* ~ u„ 


Figure 4. 2. 2. -General SPT Algorithm: Multi-Time Scales 







denote the adjoint variable of in the reduced order problem. Note that 
and Xq have the following functional dependence: 

I (2.87) x^ - x^(?Q»t^) 

(2.88) Xq = ^Q(?Q9t^) 

1 Let L« denote the cost criterion L in (2.86). Set i = l and go to the next step, 
(iii) Solve the i— initial boundary layer problem: 


(2.89) 

.dx' . 

dx 

f^(Ci_l.x'‘,Ui) 

(2.90) 

0 = 

(^^_^,x\u^) 

(2.91) 

x’(0) 


(2.92) 

x’(~) 

~ X (^^. _^jtf) 

(2.93) 

0^(u,.) 

r “ 


*th 

where the cost criterion L^. for the i— initial boundary layer problem is defined 
recursively in terms of as 




i-1, 


(iv) If i < r, denote the first component of the pseudocontrol u. calculated 

1 ^*1 ^ i 

in step (iii) by X and let X- denote the adjoint corresponding to x in the 

i— boundary layer problem. Update i by 1 and return to the beginning of 

step (iii). 

Stop if i = r; u=u^ gives the SPT approximate optimal feedback control for 

the current state 6 with time-to-go t^. 

* 

The procedure for calculating the 0— order SPT approximation of the optimal 
control law consists of solving a hierarchy of lower dimensional optimal control 
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problems. At the i = 0 level of the hierarchy, which corresponds to the slowest 
component of the state variable, we must solve the reduced order control problem. 
This is an Oq dimensional problem with finite terminal time, which may be free 
or fixed. At the subsequent levels for i > 0, which correspond to faster com- 
ponents of the state variable, we must solve the i— initial boundary layer 
control problem. This is an n^. dimensional problem with infinite terminal time, 
i.e., it is an infinite horizon control problem. 

At each level i>0 of the hierarchy the cost criterion and the asymptotic 
terminal condition for the control problem is defined recursively from the pre- 
vious level by the relation (2.94). Thus, the cost control problem at the i— 
level of the hierarchy depends on the components for 0< j< i of the current 
state 

Once the final level of the hierarchy is reached, the SPT approximation to 
the optimal' control for the current state C and the time-to-go t^ is the optimal 
control u = u^ one obtains from this level. 

4.3 Computational Aspects of the SPT Approximation 
4.3.1 Computational Difficulties of SPT Approximation 

In Subsection 4.2.3 we decomposed the original optimal control problem 

(2.1) - (2.4) of state dimension n into r+1 subproblems of dimension n. for the 

*/"h ^ ■f'h 

i~ subproblem, i = 0,l,...,r. Since we will have n^. < n for each i, the i — 
subproblem should be computationally easier to solve than the original n dim- 
ensional problem. However, the r+1 subproblems are coupled in such a way that 
the i— problem requires results from the (i-1)— problem to define the cost 
and terminal condition in the i— level problem. If we define the integers N. 
as the sums 

(3.1) N. = t n. 

’ j=0 J 

then the control subproblem at level i depends on N^.+l parameters (or N.j para- 
meters if the terminal time for the original problem is free). These para- 
meters consist of the current stat^ components for 0<j<i and the time-to- 
go t^. At the r— level, the final level, the control subproblem depends on 
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N^+l = n+l parameters— namely, the entire current state ? and the time-to-go 
t^. Thus, the final level of the control subproblem depends on as many parameters 
as the original problem and may be as difficult to solve. Moreover, it is es- 
sential to compute this final subproblem in order to obtain a realistic feed- 
back control for the original control problem. For if we stop the computation 
process at a level i<r, then we implicitly will be assuming that we can directly 
control the states for j>i. In practice this means that the u-controls 
resulting from the i-- level will not provide feedback correction for these 
higher level states for j>i. 

Thus, the SPT approximation presents us with a dilemma. To obtain a feed- 
back control law that will correct deviations from optimal of all state com- 
ponents we need to compute the final level of all the subproblems up to and in- 
cluding the final level. But to obtain a computational advantage over solving 
the original problem exactly (without SPT) by using the SPT approximation we 
need to stop our computation at some level before the final level. In the fol- 
lowing three subsections we discuss three different strategies for making the 
SPT algorithm computationally efficient. 

4.3.2 Complete Time Scale Separation 

One possibility for making the SPT algorithm more efficient is to subdivide 
the original problem into a large number of small dimensional subproblems and 
to develop fast algorithms for solving these small dimensional optimal control 
problems with general cost criteria and general terminal conditions. The 
extreme case of this procedure is to choose r = n-l so that = 1 for each 
i = 0,l,...,r. This is the approach Calise (1977, 1978) has taken in construc- 
ting feedback controls for aircraft by means of singular perturbation theory 
methods. Making each subproblem one-dimensional allows one to substitute a 
sequence of static optimization problems for the original dynamic optimization 
problem. The drawback of this procedure is that the assumption that each one- 
dimensional component of the state operates on its own time scale will be un- 
realistic when some components of the state are closely coupled together. 

The transformation of the dynamic control problem to a static optimization 
problem is based on the following simple observation. Suppose that we have a 
dynamical system with a one-dimensional state variable x and dynamic equation 
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given by 

(3.2) ^ = f(x,u) 


and cost criterion given by 


(3.3) J( 


u) =r L{x,u] 


The minimum principle asserts that the optimal adjoint variable \* for a given 
state X is chosen so that we have 

(3.4) c = min {L(x,u) + AJf(x,u)> 
u ^ 

where c is a constant not depending on x. From (3.4) we see that if u is a 
control such that f(x,u)>0 then 


(3.5) 


c - L(x,j/ 


and if u is a control such that f(x,u)<0 then 


(3.6) 


c - L(x,u] 


Defining m(x) and M(x) so that for each x we have 

(3.7) m(x) = max • f(x,u)>0^ 

(3.8) M(x) = : f(x,u)<oj 


we obtain the following upper and lower bounds on the adjoint X*: 
(3.9) m(x) ^ X* ^ M(x) 


If the optimal control u* for the given state x is such that f(x,u*)seO, 
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then the adjoint X* must be equal to either its upper or lower bound in (3.9) 
and the optimal control u* is found by solving the corresponding minimization 
(3.8) or maximization (3*7). However, for some states x it may be optimal to 
choose u* such that f(x,u*) =0, and in that case the minimum principle (3.4) 
gives no information on the adjoint variable X*. Nevertheless, the condition 

^ ■ ■ I 

that the state x tends asymptotically to a terminal state d implies that 
f(x,u*) = 0 for optimal u* if and only if x = to. Note that such an asymptotic 
terminal condition is assumed for each initial boundary layer problem in the 
SPT approximation. Note also that without assuming such an asymptotic terminal 
condition it may be possible to achieve a lower cost (3.3) with a state trajec-' 
tory which does not tend to any limit, e.g., a periodic trajectory. 

Therefore, let us add to the control problem (3.2) and (3.3) the asymptotic 
terminal condition 

(3.10) l.im x(t) = 0 ) 

Then if u*(x) denotes the optimal feedback control for state x, we have 
f(x,u*(x)) = 0 if and only if x = 6o. In this case, u*(x) is found from solving 
the static problem (3.7) when x < 03 and from solving the static problem (3.8) 
when X > oj. The corresponding optimal adjoint is equal to m(x) if x < w and Is 
equal to M(x) if x>w. This procedure for determining the optimal u and in 

terms of x is illustrated in Figure 4.3.1. 

To see how the preceding considerations apply to the multi-time scale 
control problem such as we have illustrated in Figure 4.2.2, let us assume that 
each of the initial boundary layer control problems in Figure 4.2.2 is one- 
dimensional, i.e., assume that n^ = 1 for each i = l,2,...,r. Then the initial 
boundary layer problems can be solved as a sequence of static optimization prob- 
lems as illustrated in Figure 4.3.2. In Appendix 4.2 of this chapter we show 
that this static optimization problem has a solution under very general conditions. 

Although we avoid having to solve any TPBVP's or dynamic trajectory optimi- 
zation problems fo-^ the initial boundary layer problems, the complete separa- 
tion of time scales in the boundary layer problems may result in a static op- 
timization problem which is no less difficult to solve than the corresponding 
dynamic optimization problem. In particular, note that as before the cost 
criterion and the terminal condition at each level depend on the computations of 

87 


i 


I III III 


I II I Mil I II nil I 



Figure 4.3. 1 .-Feedback Control Algorithm for One-Dimensional State. 
















the previous level. Thus, the cost criterion and the terminal condition as 
functions of become more and more complex as i increases. In the highest 
level, i = r, we will have a complex nonlinear prograrnning problem whose cri- 
terion function depends on the current state, and cannot be determined 
until the previous r-1 initial boundary layer problems have been solved. 

4.3.3 Suboptimal Approximation of the Reduced Order Solution 

A second possibility for increasing the efficiency of the SPT algorithm 
is the use of suboptimal approximations for the reduced order solution. As 
shown in Figure 4.2.2, the reduced order control problem is the first step in 
the SPT algorithm for approximating the optimal control law. By approximating 
the solution of the reduced order problem, we may be able to solve a reasonably 
large dimensional reduced order problem in an efficient manner arid leave a small 
dimensional initial boundary layer problem to solve in the next step of the SPT 
algorithm. In this way, we may reduce the computation time for the whole al- 
gorithm and at the same time, because we can consider a larger dimensional 
reduced order problem rather than an artificially small dimensional reduced 
order problem, we may also improve the accuracy of the SPT approximation. 

For simplicity let us only consider the suboptimal approximation of the 
reduced order solution in the case of two time scales as in Subsection 4.2.2. 
There is no difficulty in applying the same procedure to the multi-time scale 
case. In order to go from the reduced order control problem to the initial 
boundary layer problem (see Figure 4.2.1), it is necessary to have the adjoint 
X* of the reduced order state variable x. Not only does this adjoint define the 
cost criterion for the initial boundary layer problem, it also determines the 
asymptotic terminal condition y* of the initial boundary layer problem as a 
solution of the minimization: 

(3.11) min {L(C.y,u) + X*f(?,y,u): g(^,y,u) = 0} 

Unfortunately, we may not have the adjoint A* immediately available. For ex- 
ample, we might be given a suboptimal feedback control u^(x) and y^{x) derived 
from heuristic considerations or from previous experience and have no corres- 
ponding adjoint. Moreover, if the reduced order control law is not optimal, 
it is not clear what the adjoint of the suboptimal reduced order state should be. 


90 



To resolve these difficulties we must consider the behavior of the initial 
boundary layer solution. As we noted previously, this problem is an infinite ' 
horizon control problem with the infinite time cost criterion defined in (2.74). 
Assuming that the boundary layer state y has an asymptotic limit as in (2.73), 
then it is not hard to see that this limit must be the first component y* of 
the solution of the minimization problem in (3.11). That is, y* is the optimum 
steady state control for the reduced order problem with cost (2.74). Thus, we 
should choose the adjoint X* so that the solution of (3.11) gives a reasonable 
suboptimal control for the reduced order control problem. 

For example, if we are given the suboptimal control u^,y^ for the reduced 
order problem and we wish to find an adjoint which gives back this particular 
suboptimal control , then we must find A* such that Ug,yg is the solution of 
(3.11). For example, if u^ and y^ do not occur on a constraint boundary of y 
or u, then we may determine the adjoint X* by solving the simultaneous linear 
equations 

^ ^ ^x ^ ^ ^y ^ ^ 

( 3 . 13 ) (S-Ys-Uj) + X* (s-y^.u^) + X* fa (5^3,53) = 0 

Note that unless the dimension of the x state is equal to the dimension of the 
u control , there will not be the same number of equations as unknown adjoint 
variables in (3.12) and (3.13). Thus, in general we must expect that the sim- 
ultaneous equations in (3.12) and (3.13) will have either no solutions or in- 
finitely many solutions. If the equations (3.12), (3.13) have no solution, 
then we cannot obtain X* with this method and we must find some other way of 
choosing an appropriate adjoint X* for the reduced orddr control problem- 

Instead of trying to find an adjoint \* which will yield a given reduced 
order control y^.u^ as the solution of (3.11), let us find an adjoint which will 
yield a better reduced order control than the control ^ 5 *^ 5 • It is clear that 
such a method would be preferable to the first method described above, but it 
is not clear how we can find the desired adjoint. Fortunately, Bellman's (1954, 
1957, 1961) technique of monotone approximation provides exactly the solution we 
are seeking. 

Since monotone approximation is an application of the ideas of dynamic 
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programming, we start by briefly describing the role of dynamic programming in 
optimal control. The reader should refer to Bellman (1957) for details. Sup- 
pose that we have an autonomous dynamic system 

(3-14) ^=f(x.u) 

with the cost criterion 

(3.15) J(u) =r L(x,u) dt + C(x(T)) 

-^0 


and a fixed terminal time T. Let the state have the following initial condition 

(3.16) ., x(0) = ^ 

Note that for mathematical convenience we have fixed the terminal time T and we 
have omitted a terminal condition in favor of the terminal cost term C in the 
cost criterion (3.15). 

Let V*(^,T) denote the minimum possible cost (3.15) over all admissible 
controls u. That is, let 

(3.17) V*(^,T) = min J(u) 

u 


We can interpret V*(^,T) as the minimum cost-to-go from the state ^ with time- 
to-go T. The minimum cost-to-go V* satisfies the following functional equation: 


(3.18) Iy* = min {L(C,u) + ||* f(C«u)} 


with the boundary condition 

(3.19) V*(C.O) = C(C) 

If u*(^,T) is an optimal control to use at state ^ with time-to-go T, then (3.18) 
impl ies 

(3.20) 1^* = L(5.u*(5,T)) + II* f(5.u*(S,T)) 
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Monotone approximation uses the boundary condition (3.19) and the linear equa- 
tion (3.20) in V* to find a cost-to-go function corresponding to a given 
suboptimal control law u^. Thus, suppose that satisfies 

(3.21) ^ = L(^.U3(^.T)) + ^ f(^,Ug(^,T)) 

with the boundary condition (3.19). Then if we generate a new suboptimal control 
Ug by minimizing the function 

3V 

(3.22) L(?,u) + ^ f(5,u) 

with respect to u, the new cost associated with this new control law u^ will be 
less than the cost for the control law u^. By repeating this process one obtains 
a sequence of cost-to-go functions which decrease monotonical ly to the optimal 
cost-to-go V*. 

To obtain an appropriate adjoint to use in the initial boundary layer prob- 
lem, we only need to apply the monotone approximation method once. Thus, given 
a suboptimal control law u^, we solve (3.21), (3.19) to obtain Vg(^,T). The 

adjoint corresponding to the current state ^ with time-to-go T is then the partial 
8V 

derivative . Thus, we choose as an adjoint the variation with respect to 
the initial state or the cost-to-go from that state. 

To summarize, the monotone approximate method for obtaining the reduced 
order adjoint consists of the following steps: 

(i) Solve the linear partial differential equation 

(3.23) 1^ - 

with an "appropriate" terminal condition of the form 

(3.24) V(C,0) = C(?) 

If a terminal condition x^ is specified for the reduced order control problem 
rather than a terminal cost C as above, then we may introduce a fictitious cost 
such that C(x) =0 for x=x.p and C(x)=+°°. A finite cost in (3.24) would be more 


93 



tractable mathematically and computationally, and perhaps more realistic than 
a fixed terminal condition. After solving (3.23), (3.24), obtain the adjoint A* 
from 

(3.25) ys.T) = II (5,T) 

(ii) Using the adjoint calculated in step (i), minimize the function 

(3.26) L(C,y,u) + A^(C,T)f(^,y,u) 

with respect to the reduced order controls 0 and y. The solution of this mini- 
mization problem provides the necessary asymptotic terminal condition for the 
initial boundary layer control problem. 

To. implement this method for finding the reduced order adjoint, we must ob- 
tain an efficient solution to the partial differential equation in (3.23). 

Solving partial differential equations is never an easy task, and the computa- 
tion of the adjoint A^ would most likely have to be done off-line. Nevertheless, 
note that for a given suboptimal control the solution of (3.23) will be much 
easier than the solution of the nonlinear dynamic programming partial differential 
equation (3.18) for the optimal cost-to-go. Moreover, a number of techniques 
are available for solving and approximating the linear problem. 

In Figure 4.3.3 we have illustrated the basic procedure of using a sub- 
optimal reduced order control. Note that in the multi-time scale case we may 
use this suboptimal strategy at any level of the SPT algorithm. Thus, we may 

+ U 

compute a suboptimal solution of the i— initial boundary layer control problem 
and use one of the above methods to compute an adjoint variable to use in de- 

‘f'h 

fining the cost criterion for the (i+1)— problem. 

4.3.4 Linearization of the Boundary Layer Problems 

A third possibility for increasing the efficiency of the SPT algorithm is 
the linearization of the initial boundary layer control problem around the 
nominal provided by the first component of the pseudocontrol of the reduced 
order solution. Implicitly, this approximation assumes that the "fast" boun- 
dary layer variables are approximately in equilibrium; that is, equation (2.54) 
is approximately correct. As we have done in the previous two subsections, we 
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Figure 4.3.3, -SPT Algorithm Using Suboptimal Reduced Order Solution. 









will first discuss the two time scale case in detail and then briefly discuss 
the multi-time case. 

Consider the two time scale problem in Subsection 4.2.2 and let the current 
X state be fixed at Let X*, y*, and u* be the corresponding optimal adjoint 

y\ 

X^ and pseudocontrols y, u to use for this value of the x state in the reduced 
order control problem. We are going to linearize the boundary layer control 
problem (2.71)- (2.74) around these optimal reduced order values. In the 
following derivation let us denote quantities evaluated at the optimal reduced 
order solution by a bar, e.g., Let 6y denote y - y* and let 6u denote ti - u*. 

ay 

Expanding (2.71) to first order gives 

(3.27) 6y’=||6y+||5u 

where the derivative with respect to x is denoted by a prime, i.e. , ( )' = - • 

Let us define matrices A and B as 

(3.28) A = |a (5,y*,u*) 

(3.29) B = |a (5,y*,ii*) 

Then (3.27) represents the time invariant linear system 

(3.30) 6y' = A6y + B<Su 

To obtain the quadratic cost criterion for the linearized problem we first 
define the Hamiltonian for the boundary layer problem, that is, 

(3.31) H(y,u) = L(^,y,u) + A*f(C>y>u) + X^g{^,y,u) 

From Bryson and Ho (1975) we see that the cost criterion for the linearized 
boundary layer problem is defined in terms of the following three matrices which 
consist of the second order derivatives of this Hamiltonian. 

(3.32) Q = 

9y 


96 



(3.33) 


2 

= ^ H 
3y3u 

(3.34) S = ^ 
3u 


In the expressions (3.32) - (3.34) the bar denotes that the expression is evalua- 
ted at the reduced order nominal value. Thus, for example, we have 


(3.35) 


^ ^ (e.y.u) (5.y.u) + L ^ (s.y.u) 

3y 3y * 3y^ ^ 3y 


Evidently, we need to calculate the Lagrange multiplier in the reduced order 
problem if we are going to linearize. For example, we can accomplish this by 
solving the linear equations 


( 3 - 36 ) 


Note that some of these equations will be redundant for the purpose of obtaining 



Using the expressions (3.32) - (3.34), we can write the cost criterion for 
the linearized problem as 


(3.38) 


r 

Jo 


h [5y 


Q 6y + 


26y^ R 6u 


6u^ S 6u] dx 


Hence, the linearized approximation of the boundary layer control problem gives 
a time invariant, linear quadratic regulator problem. The feedback control is 
simply calculated by solving a system of quadratic equations (the so-called al- 
gebraic Riccati equation) to obtain a constant gain matrix G. For this problem 
G is obtained by solving the following equations for G and K simultaneously: 

(3.39) G = -[S‘3(KB+R)]3^ 
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(3.40) Q + 2RG + G^SG + KA + KBG = 0 

The optimal control law for 6u in terms of 6y is then 

(3.41) 6u = G6y 

To summarize, the procedure for finding the linearized boundary layer control is 
to compute the matrices A, B, Q, R, S from (3.28), (3.29), (3.32), (3.33), (3.34) 
for the current value of the reduced order state, and then to solve the equations 
(3.39), (3.40) to obtain the gain G which gives the boundary layer control law 
in (3.41). 

In the multi-time scale case we may carry out the linearization at any 
level of the SPT algorithm. Thus, it may be appropriate in some cases to lin- 
earize the control problem for the i^ control subproblem around the nominal 
values proved by the (i-1)— control subproblem. Other variations are also 
possible. For example, we may decide to linearize the boundary layer problem 
described above around y* only and maintain the nonlinear dependence of the boun- 
dary layer problem on the control u. In any case, the decision whether or not 
to linearize is determined by the "error" between the actuaT value of a variable 
on the i^ level of the SPT algorithm and its "optimal" value on the (i-1)^ 
level. By monitoring this error an algorithm can decide when to switch from a 
more accurate and also more difficult nonlinear control to a more, efficient 
linear approximation. One such algorithm is illustrated in Figure 4.3.4. 

Before concluding this section, let us note some of the advantages and dis- 
advantages of the linearization procedure described in this subsettion. The 
major disadvantage 'is that we are usi-ng a linear model for a nonl'ineor system. 

In particular, this approximation will be valid only when the boundary layer 
states are sufficiently close to their "optimal" values calculated in the re- 
duced order solution. The major advantages of linearization are that (1) the 
linear control is easier to compute than a nonlinear control, and (2) because 
the linear control is easier to compute, we may treat a larger dimensional boun- 
dary layer control problem than otherwise: thus, we may avoid having to intro- 

duce multiple time scales artificially in order to obtain a tractable computation. 
Note that once a linear problem is obtained it is possible to apply the General- 
ized Multiple Time Scales method (GMS) of Ramnath and Sandri (1969) to improve 
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the efficiency of computation. The time scales for the linearized problem are 
determined by checking the eigenvalues of the linear system and thus, we can 
naturally determine the perturbation parameters. This multiple time scaling of 
the linear system decomposes the control problem and reduces the computation 
required. In addition, the time scaling improves the numerical accuracy of 
computation. This latter aspect is most important in view of the numerical sen- 
sitivity of the linearization. This sensitivity arises from the presence of 
tabular, non-analytical functions (aerodynamic coefficients) whose derivatives 
must be given for the linearization. In Figure 4.3.5 we have summarized the 
principal advantages and disadvantages of the three strategies we have discussed 
in this section for making the SPT algorithm more efficient. 


4.4 Accuracy of the SPT Approximation 

4.4.1 State Space Dependence of the Accuracy of the SPT Approximation 

The SPT algorithm discussed in the previous sections gives a control law 
u(C,T) for the control problem (2.1) - (2.4) in terms of the current state ^ and 
the time-to-go T=t^-tQ. Let J(u;^,T) denote the cost-to-go starting at the 
state ^ with time-to-go T and using the control law u. If V*(^,T) is the mini- 
mum cost-to-go from the state ^ with time-to-go T, then the expression e(u;C»T) 
given by 

(4.1) e(u;5.T) ^ J{u;5,T) - V*(5,T) 

defines a measure of the accuracy of the control law u compared to the optimal 
control. Since V* is the optimal cost, it is clear that e is always nonnegative. 
However, as the initial state ^ and the time-to-go T vary, we may expect that 
the error e will vary also, being greater in some regions and smaller in others. 

Let u specifically denote the SPT approximate optimal control law. For 
any fixed initial state ^ ^ w and time-to-go T the error e(u;^,T) depends on the 
singular perturbation parameter e. As e decreases to 0 we expect the error to 
approach 0. However, the error does not approach 0 uniformly in i and T. That 
is, for a fixed perturbation parameter e, the error may vary greatly for differ- 
ent values of ^ and T. In this section we would like to show that the SPT ap~ 
proximatlon error Increases as the Initial state ^ approaches the terminal state 
0 ). That Is, as one nears the terminal target state, the SPT approximation breaks 
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SUMMARY OF PRINCIPAL ADVANTAGES AND DISADVANTAGES OF 
COMPLETE TIME SCALE SEPARATION, SUBOPTIMAL REDUCED ORDER SOLUTION, 
AND LINEARIZATION OF BOUNDARY LAYER CONTROL PROBLEM 


1. Complete time scale separation (see Subsection 4.3.2) 

Advantage : Reduces dynamic trajectory optimization 

problem to static optimization problem 
for control law. 

D isadvantage : Introduces multiple time scales arti- 

ficially and will be inaccurate if some 
state variables are highly coupled; the 
sequence of static optimization problems 
may be as difficult to solve as the ori- 
ginal dynamic problems for a large dimen- 
sional problem where many time scales are 
necessary. 

2. Suboptimal Reduced Order Solution (see Subsection 4.3.3) 

Advantage : Allows one to avoid solving the reduced 

order problem exactly; thus, it is pos- 
sible to treat larger dimension reduced 
order problems. 

Disadvantage : The suboptimal adjoint function must 

be computed and approximated off-line, 
although this should be easier to do 
than to compute the optimal adjoint. 

3. Linearization of the Boundary Layer Problem (see Subsection 4.3.4) 

Advantage : Linearized controls can be computed ef- 

ficiently and larger dimension boundary 
layer problems can be treated. 

Disadvantage : Linearization is valid only when the 

boundary layer state is near its optimal 
reduced order value; thus, linearization 
essentially gives only a control to track 
the reduced order solution, and it will 
be able to do this only as long as the 
pseudocontrol for the reduced order prob- 
lem does not jump discontinuously. 

Figure 4.3. 5. -Comparison of Advantages and Disadvantages of Complete Time 
Scale Separation, Suboptimal Reduced Order Solution and 
Linearized Boundary Layer Approximation. 
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down. Before discussing this breakdown in more detail, however, let us note that 
the opposite phenomenon is true when the initial state is far from the terminal 
target, at least in some cases; that is, the SPT app2*oootmatton 'improves as the 
init'lai state ^ increases its distance from the terminaZ state o). 

As we will see in the next chapter, the aircraft dynamical system has the 
general form 

(4.2) ^=f(y,u) 

(4.3) ^ = g(y,u) 

where the right hand side of (4.2) and (4.3) do not depend on x. In the aircraft 
problem, (4.1) corresponds to the equations for the horizontal position coordi- 
nates of the aircraft. The variables y are such that either through direct 
state constraints or by definition (as in the case of angles) each component of 
y has a maximum value. In the aircraft example y includes height, velocity, 
heading angle and so on. The x variables, on the other hand, have a potentially 
infinite range. If we scale the variables x and y such that the scaled quan- 
tities x‘ and y' attain a maximum variation of order of magnitude 1, and if we 
scale the time variable t so that the scaled time t‘ varies from 0 to 1, then 
the original dynamic system (4.2), (4.3) is transformed to 

(4.4) ^=f'(y',u) 

(4.5) = g'(y' ,u) 

where the parameter e is essentially the ratio of the maximum variation of y to 
the initial distance from x to its target state x.p. Moreover, the functions f 
and g' in the aircraft example do not depend on the parameter e. Thus, in the 
minimum time problem for the aircraft eocample there is at least one natural 
singular perturbation parameter, and this parameter is proportional to the 
distance from the target. Since we expect that the SPT approximation is best 
when £ is small, we expect that the SPT approximation for the aircraft will be 
good when the range is sufficiently long, but that the SPT approximation will 
break down near the target when the range is short. Note that for the example 
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we have considered, the SPT approximation starts to break down when the scale 
parameter e in (4.4), (4,5) increases above ,07. This corresponds to target 
distances less than 163 to 176 Km. 


4.4.2 Inaccuracy of the SPT Approximation Near the Terminal Target State 
It is the purpose of this subsection to show that one must be especially 
careful when using SPT approximations to obtain closed loop, control law 
approximations rather than open loop, control trajectory approximations. At the 
terminal time the SPT control law will not be the asymptotic approximation of 
the actual optimal feedback control law. To illustrate what goes wrong we will 
consider the following simple linear quadratic control problem. 

(4.6) ef = u 

with the initial and terminal conditions x(0)=^ and x(T) = 1, respectively . The 
cost cr iter ton for the problem is 


(4.7) J(u) = + dt 

-’o 

The exact optimal control law is easily found to be 


(4.8) 


’ (l_e-2T/e) 


The 0^ order SPT approximation, the first term of (2.6), is also easily found 
to be 


(4.9) u(e,T) = -x+e“^/^ 

Note that we have written "u" in (4.9) to indicate that both the u and u (as 
well as u) terms are accounted for. 

To see how the SPT approximate control law in (4.9) differs from the exact 

optimal feedback control law given in (4.8), we approximate the latter to first 

order. First consider the region in which the time-to-go T satisfies T»c. 

-T/c 

The optimal feedback law f4.8) is then approximately, to order 0(e ), 
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(4.12) Q(C,T) - -(C-1) 

In this control problem (4.6), (4.7) one desires to reach the terminal 
target x=l by a fixed terminal time, say t^=l. Then as the current time t 
nears the terminal time t^ the optimal control law (4.8) takes this into account 
through the time-to-go T-t^-t becoming small. As the time-to-go T becomes 
smaller, the optimal control law places more emphasis on hitting the target 
exactly at the right final time t^ and less importance on reducing the quadratic 
cost. The SPT control law (4.9), on the other hand, resembles a linear regula- 
tor control around the terminal state x = l. Although such a control would steer 
the state toward x = l, it would not be able to do so in a finite amount of time. 
Hence, at the terminal time t^ the SPT controlled state will miss the desired 
terminal state. To see more clearly what trajectory the SPT control produces, 
let us integrate (4.6) with u = u. Then we find that if C is the initial state 
with time-to-go T, the final state x(t^) using SPT control will be 

(4.13) x(tf) = + I 

If x = 1 is the desired terminal state, then the. error between this and the actual 
terminal state (4.13) is given by 

(4.14) 1 - x(t.p) = 1 - f + I 

Thus, 1f the initial state C is such that is small and the initial time- 

to-go T is large compared to e, then the final error at the terminal time t^ 
will be approximately 

(4.15) 1 - x(t.p) = 1 - f 

In fact, for any initial condition if the initial time-to-go is large enough, 
then (4.15) will be true. 

The situation is worse than indicated by the approximation of u given in 
(4.12)! If we use the SPT control all along and the state x approaches too 
closely to its optimal reduced order value x*=0, then x will still be close to 
the reduced order value at the final time. To be exact, if at any time we ever 
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have X < then x(t^) < e. 

A similar phenomenon occurs for more complex systems than the simple one 
we have examined here. Essentially, the SPT control follows the reduced order 
solution and is unable to anticipate discontinuities in the reduced order solu- 
tion in order to cross the jump of the discontinuity in the required amount of 
time. In this example the reduced order value of x was 0 and the desired ter- 
minal value was 1. Although the SPT control law begins to steer x away from 
the reduced order value 0, it does not do it fast enough to meet the terminal 
value and the resulting error (4.15) is fairly large. 

The trouble we have described is not limited to the terminal region: so- 

called internal boundary layers exhibit a similar behavior. To describe the 
situation, let us refer to the two time scale problem of Subsection 4.2.2. In- 
ternal boundary layers occur when the reduced order control problem gives a 
pseudocontrol y*(?,T) with discontinuities for some values of current x-state ^ 
and time-to-go T. Such discontinuities indicate that the real y-state, which 
must be continuous, changes very rapidly. This rapid change is approximated 
by a boundary layer and this requires us to modify our basic time scale approxi- 
mation (2.38) of the optimal control u* for the full system as follows: 

00 

(4.16) u* ^ YL {[U|^(t) - Uj^(tQ)] + [U|^(p) - Uj^(t^)] + Uj^(t) 

0 

+ [U|^(a) - U|^(t^)]} 

where t^- is the time at which the reduced order solution exhibits a discontinuity 
in the pseudocontrol y* and p is the fast time scale for the internal boundary 

t-ti 

layer given by — • 

The argument we made in Subsection 4.2.1 concerning the terminal boundary 
layer also holds for the internal boundary layer — as long as the time t is not 
close to the time t^. . We may neglect the contributions Uj^ from the internal 
boundary layer. Near to the internal boundary layer time the SPT approximation 
breaks down in the same way that it did at the terminal time. Essentially, the 
SPT control law does not anticipate the sudden rapid change in the y variable. 

An optimal control would begin steering y across the discontinuity before the 
internal boundary layer time t^. was reached. The SPT control, on the other 
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hand, does not do this until t^ is reached, and then it steers y to the reduced 
order value y* after t^. 

At the present time we do not have a solution to this fundamental problem 
of using SPT control law approximations near the terminal target or at internal 
boundary layers, but clearly the problem is essential for any application of 
SPT methods to obtain feedback control laws. At the very least, it will be 
necessary to identify the regions where the SPT feedback law breaks down so that 
one may switch from the SPT control to a better control law. For obtaining the 
control law near the terminal target, it may be possible to develop useful ap- 
proximations different from the SPT approximation. For example, linearizing 
some of the state variables around their terminal values may prove useful. How- 
ever, whatever method of approximation is used near the terminal target, the 
important problem to solve will be to determine how to match this approximation 
with the SPT approximation — in other words i to determine when and how to switch 
from one type of control to another. 
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APPENDIX 4.1 

MULTI-TIME SCALE SPT FORMULATION AND CALCULATION 

In this appendix we present the SPT formulation and calculation of the 
feedback law in more detail than in Subsection 4.2.3. Throughout this appendix 
we will maintain the same notation as we did in that subsection. In addition, 
let us define as 

(1.1) x^. = (x^,...,x’’) 

just as we defined in (4.2.79). Likewise, define f^. by 

(1.2) f^Cx,u) = (f^(x,u),. .. ,f^(x,u)) 

The SPT method seeks to approximate the control system in equations (4.2.1) 
through (4.2.4) by finding an asymptotic expansion to the singularly perturbed 
system 

(1.3) e. = f’(x,u) 
with the same cost criterion 

(1.4) J(u) =r ^ L(u,x) dt 

-^0 

the same initial conditions 

(1.5) x’(tg) = e’ 

and the same terminal conditions 

(1.6) x’(tf) = J 
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In (1.3) the parameters e. are such that Eq- 1 and for each i we have 

(1.7) - 0 

as £^.->-0. In reality, of course, we will be using the approximation for 
for all i. The relation (1.7) represents the assumption that the variables 
vary on different time scales which become increasingly faster as i increases. 

Although the multi -time scale SPT approximation is conceptually no more 
difficult than the two time scale case, the notational difficulties in the multi 
time case become extreme. Therefore, we will try to simplify the problem at 
an abstract level as much as possible before working on the control problem. 
Working at the abstract, more general level we can avoid some of the notational 
obscurities and maintain a modest level of clarity. 

The essential idea of the SPT approximation is to approximate a trajectory 
z(t) for tQ<t<t^, where z represents all the state components, adjoint com- 
ponents and controls. The 0^ order approximation has the general form: 

(1.8) z(t) =: [z(t) - z(tQ)] + z(t) + [z(a) - z(tf)] 

where the vectors t and a represent the multi -time scales, namely, 

(1.9) T = (Tq,Tj,...,T^) 

(1. 10) o = (aQ,a^,. . . ,a^) 

t-t^ t-t^ 

where t. is defined t . = ^ and a. is defined a . = As for the two time 

1 1 £• 1 1 e. 

scale case, we need only consider the initial boundary layer represented by the 
term z in (1.8). For multiple time scales, the 0^ order approximation of the 
initial boundary layer term has the form 

(1.11) z(?)=i^(xj)+2[$j,^(xj,j)-zj(0)] 

where the terms z.j do not denote different components or combinations of com- 
ponents of z as we defined x^ and x^. previously. The z- all have the same 
dimension and represent the approximations of z on each of the different time 
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scales T^.. 

The expressions z^- must have the asymptotic property that 

(1.12) T-IjS = 2i_i(0) 

for each i>l. For this reason, the approximation (1.11) of z(0) is given by 
(1.12) z(0) - yo) 

Compare this to our discussion in Subsection 4.2.1. The approximation (1.11) is 
analogous to the first term of (4.2.44). Thus, to obtain the SPT approximation 
to z(0) we need to calculate the value of z^(0), the fastest initial boundary 
layer approximation at t^=0. To do this, however, requires first computing 
the slower initial boundary layer approximations. 

The dynamic system for the time scale t. is obtained by transforming (1.3) 
to the independent variable t^. Thus, we have 

fly j ^ -j 

(1.14) I7 = / f^(x,u), Os j<1 

i i 

(1.15) = f’(x.u) 

i 

(1.16) = f’(x,u), raj>i 
i 

4.^ 

The 0— order approximation is obtained by letting and using the relations 

(1.7). The equations (1. 14) - (1. 16) then become 

(1.17) = 0 

i 

(1.18) = f’(x,u) 

^ (1.19) 0 = H(x,u), raj>1 
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The equations (1.20) and (1.21) represent the dynamical equation and equality 
constraint, respectively, for the i— level system operating on the time scale 
T^.. The initial conditions for the state in this problem are those given 
by (1.5). The terminal conditions, however, are determined as a result of the 
solution of the control problem for the (i-1)— level. These terminal conditions 
are chosen to satisfy the asymptotic relation (1.12). Likewise, the cost cri- 
terion for the i— level is determined on the (i-1)^ level by the adjoint for 
the (i-1)^ level. The terminal condition is given in (4.2.92) and the cost 
criterion is given in (4.2.93) of Subsection 4.2.3. 
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APPENDIX 4.2 

THE EXISTENCE OF SOLUTIONS TO THE BOUNDARY LAYER CALCULATION 
IN THE CASE OF COMPLETE TIME SCALE SEPARATION 

In this appendix we consider the SPT algorithm's calculation of the control 
law for a dynamic system for which we have assumed complete time scale separation. 
To simpl ify matters let us first consider the case of two time scales where x 
and y represent the slow and fast variables respectively. We assume that both 
are scalar state variables. Let u denote the vector control variable for the 
problem, and let L(x,y,u) denote the integrand of the cost criterion. The dynamic 
equations for the system are 

(2.1) ^=f(x,y,u) 

(2.2) ^ = g(x,y,u) 


As discussed in Subsection 4.3.2, the SPT algorithm calculates the control law 
for this system in two stages as follows. Let x^ denote the given final value 
for the X state variable. The first stage of the SPT calculation is to solve 

(2.3) * g(x,y,u) = o, f(x,y,u) > 0^ = 
if X < x^, and to solve 

(2.4) : g(x,y.u) = 0. f(x,y,u) < 0 J = -A*(x) 

if X > x.p. Let y*(x) denote the optimal pseudocontrol and let u*(x) denote the 
optimal control which solve (2.3) and (2.4). 

The second stage of the SPT algorithm is to solve 


(2.5) 


min 


L(x,y,u) + X*(x)f(x,y,u) 
g(x,y,u) 


g(x,y.u) > 0 } = -X*(x,y) 


if y<y*(x) and to solve 
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( 2 . 6 ) 


: g(x,y,u) < 0 


max 


t L(x,y,u) + A*(x)f(x,y,u; 
~~ g(x,y^) 


= -A*(x,y) 


if y>y*(x). Let u*(x,y) (denote the optimal control which solves (2.5) and (2.6). 

The problem is whether or not we can actually obtain a solution to this 
algorithm, at least in principle. The question arises because it is not at all 
clear that the minimum or the maximum exist in (2.3), (2.4), (2.5), (2.6). In- 
deed, given numerical inaccuracies it is possible that numerical optimization 
algorithms will not converge to a solution unless care is taken in setting up 
the problem. It is the purpose of this appendix to show that under reasonably 
general circumstances it is possible, at least theoretically, to obtain finite 
solutions to the minimization and maximization problems above and that we can 
obtain at least one solution u*(x,y). Thus, nonconvergence of optimization sub- 
routines within the SPT algorithm must be due to numerical difficulties (which 
can be repaired) and not due to the lack of a solution to the minimization or 
maximization problems within the SPT algorithm. Thus, the SPT algorithm will 
yield a control value for any input state value--at least in theory. This fact 
was reassuring to know when some stages of the algorithm proved to be sensitive 
to numerical error. 

The first step in our argument is to prove a small lemma. Notice that in 
this lemma the optimization problem represents an abstracted and boiled down 
version of one of the stages of the SPT algorithm mentioned above. Thus, to 
find solutions to the SPT algorithm above, we will only have to apply the lemma 
to each stage, one after the other. The reader who is unfamiliar with the mathe- 
matical analysis required in the statement and proof of this lemma may skip to 
the theorem and remarks following the theorem. The analysis required here may 
be found in any introductory book on real analysis such as Rudin (1964). The 
main result of mathematical analysis which we use is the fact that a continuous 
function for a compact set K always achieves a finite minimum f(Z|) and a finite 
maximum f(z 2 ) for some z^ and Z 2 in K. 


Lemma 

Let K denote a compact (closed and bounded) set in and let f and L denote 
continuous functions from into R. Suppose that for all z in K such that 
f(z) = 0 we have 
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(2.7) L(z) > 0 


Furthermore, suppose that there is a Zq in K such that f (Z q) > 0. 

Then there exists a z* in K such that f(z*) >0 and such that for all z in 
K such that f(z) > 0 we have 



That is, z* minimizes L(z)/f(z) for all z in K such that f(z)>0. 

Proof of the Lemma 

The proof is rather simple. Essentially we show that (2.7) implies that 
as f(z) approaches 0, the ratio L(z)/f(z) approaches +«>. Thus, we can show that 
there is a positive constant 5 such that we can restrict ourselves to z such 
that f(z) > 6 is the minimization of L(z)/f(z). Once this is done, standard 
compactness arguments show that z* exists. 

Let Z denote the set of z in K such that f(z) = 0. Note that Z is compact 
and since L is continuous, L must have a nonzero minimum on Z, That is, there 
is a constant e>0 such that 

(2.9) L(z) > e 

for all z in Z. Let Z denote the set of all z in K such that the distance 

P 

d(Z,z) from z to the set Z satisfies d(Z,z) <p. The set Z^ is relatively open 
in K, and K-Z^ is compact. Since L is continuous and since Z is compact, (2.9) 
implies that for sufficiently small p, we have 

(2.10) L(z) > I 


for all z in Z . 

P 

Let Kg denote the set of points z in K such that f(z) > 6. Note that Kg is 
compact and disjoint from Z. Since Z is also compact, we can choose p small 
enough so that Kg and Z^ are also disjoint. 

Choose 6 so that 6<f(zQ) and so that 


( 2 . 11 ) 


26 
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Next choose p so that is disjoint from and so that (2.10) is satisfied for 

z in Z . Then it is clear that for z in Z we have 

P ■ ' ■ 


, L(z 
fTZn) ^ 



Thus, the minimum of L(z)/f(z) must lie in The function L(z)/f(z) is con- 

tinuous on this compact set K-Z and hence, standard theorems imply the existence 

P 

of a z* minimizing L/f. 

Given the lemma above, it is fairly easy to show that under reasonably gen- 
eral conditions, the SPT algorithm generates a well-defined solution. In the 
following theorem we note the general conditions under which a solution can be 
obtained to the SPT algorithm. 


Theorem (Existence of Solutions to the SPT Algorithm ) 

Assume the following are true: 

1) For all X there are compact sets, J and K such that yeJyCR^ 
and u e Kjj c r"'. 

2) For all x the functions (y,u) ^ f (x,y,u) , (y,u) ^ g(x,y,u) and 
(y,u) ^ L(x,y,u) are continuous. 

3) For all x,y,u such that f(x,y,u) = 0 and g(x,y,u) = 0 we have 
L(x,y,u) > 0, if X x^. 

4) The optimal pseudocontrol y*(x) solution to the first stage of the 
'SPT algorithm is unique for each x. 

Then the SPT algorithm has a solution u*(x,y). 


Remarks 

Before sketching the proof of this theorem let us comment on the above 
assumptions (l)-(4). The first assumption in (1) above merely requires that 
we can restrict y and u to a bounded set in R"^ and R^ respectively for any 
given x. The second assumption is obvious. The third assumption (3) is natur- 
ally satisfied for many cost criteria (e.g., the minimum time and quadratic 
criteria, but not the minimum fuel criterion). The first three assumptions 
guarantee a solution to the first stage of the SPT algorithm. The fourth 
assumption is necessary to guarantee a solution to the second stage of the SPT 
algorithm. If there is more than one solution y*(x), then the fast boundary 
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layer problem in the second stage of the SPT algorithm will have more than one 
possible terminal condition. Such ambiguity would require an investigation 
going beyond the first order necessary conditions provided by the minimum prin- 
ciple to obtain the actual optimal control u*(x,y) for the fast boundary layer 
problem. 

Proof of Theorem 

The first three assumptions (1) - (3) permit direct application of the lemma 
to prove the existence of a solution y*(x), u*(x) to the first stage of the SPT 
algorithm, for x^ix^. If x x^ and y?^y*(x), then the fourth assumption (41 
implies that for all u such that g(x,y,u) = 0 we have 

(2.13) L(x,y,u) + Aj^f(x,y,u) > 0 

From (2.13) we can again apply the lemma directly to prove the existence of a 
solution u*(x,y) to the second stage of the SPT algorithm. Note that u*(x,y) 
need not be unique. /// 

For more than two time scales with complete time scale separation, the exis- 
tence results are similar. The essential assumption, besides the obvious con- 
tinuity and compactness assumptions, is the positivity of the original cost 
criterion when the slowest state variable is not at its terminal value (i.e. , 
we also assume that the generalization of assumption (3) is true), and in addi- 
tion, the optimal pseudocontrol for the next faster state variable is unique at 
each stage. For example, if there were a state variable z which was faster 
than y, then we would have to assume that the pseudocontrol value for z cal- 
culated at the second stage (y boundary layer problem) was unique. Note again 
that this assumption is necessary even to make sense of the SPT algorithm. 



(1.16) L = Ljh.E)a 

V^x = Vx('’) 

Note that in the energy state formulation the velocity V is simply shorthand for 
the function of E and h given by 

(1.18) V = /2gir-h) 

Note also that E is a specific energy measured in units of height. The constraint 
bounds and a are constants in (1.9) and (1.11) respectively. The constraint 

(1.13) is present only to guarantee that the square root in (1.18) is well- 

defined. In a model which includes both E and h as state variables this con- 
straint will not be necessary, since the dynamic equations will never allow the 
height h to exceed the energy E. However, in reduced order models which include 
only E the constraint will be necessary since h will be a free control. 

The control problem is to steer the system (1.3) - (1.8) from an initial 
state (x^,y^. ,E^ ,3^. ,h^. ,Y^) at t^ to a final state (x^,y^,E:jr,3f »h^,Y^) in minimum 
time t^'tQ. Thus, the integral cost criterion for this problem is the same as 
(3.2.6), namely 


(1.19) 


J(a,a,u 



1 dt 


Since the final time t^ is not fixed, the Hamiltonian for this problem is 
identically 0 along an optimal trajectory. Note that the Hamiltonian is given 
by 


(1.20) H = 1 + A^V cos3coSY + AyV sin^cosy + A^ 

(i + uTsina)sina 

* ^ V * V 




mg 

(L + uT|^^^sina)cosa- mg cosy 
mV 


To approximate the optimal control law for this minimum time problem we will 
follow the procedure outlined in Section 4.2 of solving a sequence of smaller 
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dimensional control problems corresponding to different time scales for the 
problem. 

At this point we will indicate how the time scale separation was chosen. 
Previous work such as Kelley (1971a, b, 1973c), Kelley and Lefton (1972a) and 
Parsons (1972), which considered three dimensional maneuvers, and hence con- 
sidered turns, indicated that the basic time scale separation consisted of the 
successively faster groups (x,y), (E,3) and (h,y). Generally, 6 was considered 
faster than E and y was considered faster than h, although Kelley (i973c) noted 
that E might be considered faster than 3 in some situations,. As we have dis-; 
cussed in Section 4.4, the accuracy of a feedback control calculated with a 
particular time scale separation depends on the region of state space in which 
the control is applied. It seems clear from examples that no one time. scale; 
separation will be accurate for all regions of state space. Indeed, near the 
terminal target all time scale separations are inaccurate, i.e. , the SPT assump- 
tion of any time scale separation is invalid near the terminal target. 

Although there is as yet no systematic theory to determine the proper time 
scale separation for a given region of state space, we present the following - 
observations for guidance. The time scale separation for an optimization prob- 
lem seems to depend on two factors: (1) the sensitivity of the immediate cost 

with respect to the current state variables and (2) the sensitivity of the long 
term cost with respect to the current state variables. It appears generally 
that the most sensitive states should be considered slower than the least sen- 
sitive states. In addition, it appears that the sensitivity of the long range 
cost is important when the target is near (short term problem) and that the sen- 
sitivity of the immediate cost is important when the target is far (long range 
problem). It is important to note that sensitivity analysis of the two kinds 
of cost, immediate and long-range, can result in different choices of time 
scales. The right choice depends on the proximity of the target. 

The sensitivity of immediate cost is determined mainly by the relative 
magnitude of the time derivatives, at least in the minimum time optimization 
problem. The most sensitive states (with respect to immediate cost) are the 
ones with the smallest derivatives and hence the slowest in terms of actual 
time variation. 

We carried out a sensitivity analysis for the aircraft model presented in 
Chapter 1 as follows. We first scaled all state variables and derivatives so 
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that the state variables and time-to-go were of order one (the maximum variation 
of each variable was one). That is, we let x = ^ where = range-to-go, 

V = X where V^ = maximum speed, h = T- where h = the cruise altitude, E = -f where 

^ ^ . "s '"s- 

cruise energy, t = -i- where t_ = minimum time-to-go, namely R./V . Since $ 
s s 

and Y are angles and already vary on the order of one, we left them unsealed. 
Similarly, we scaled lift, drag and thrust in terms of their maximum values. 

With this scaling we obtained the following equations: 


Hy ”■ 

(1.21) -yf ^ C0S3C0SY 
dt 


( 1 . 22 ) 

(1.23) 

(1.24) 

(1.25) 


^ = 
dt 


V sin3cosY 


= 5pVF. 

dt E I 


dh r w • 


dy _ 
£ — ^ — 

dt 


- 6 


COSY 
W V 


d6 . 


(1.26) = = 


dt V sinY 


where F|| is the force acting parallel to the plane, namely 

(1.27) ^j| ” coset - 6|^D 

in terms of scaled thrust f and drag D. Similarly, is the perpendicular lift 
force 

(1.28) F^ = L + 6yT si not 

in terms of the scaled lift L. For the aircraft used, the parameters had the 
approximate values 6^= 1.1854, 6^ =.0430, 6p=.0313, 6y=.2503, 6j=.0475 and 
e = 3.0516/N where N = number of kilometers in the range-to-go. 
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From this scaling in equations (1.21) - (1. 26) we clearly see the separation 

of x,y from E,h,y,6 for long-range (large N) problems. In addition, note that 

. dF 

there IS a milder separation of E,h,y and 3. The order of ~ is <5^6^- .04, the 

order of -™ is 6.,^ .25 and the order of both ^ and ~ is of order 1. On the 
dt V clt dt 

basis of this scaling, there appears to be a time scale separation making x,y 
the slowest variables with E next, h after E and with y,3 together as the fastest. 
However, two important points are in order. First, in our model aircraft the 
angle-of-attack a was restricted so that a^O. Thus, in equation (1,25) we can 
only have the dominant term F^ nonnegative. Practically speaking, this means 
that the aircraft can accelerate upward much faster than it can accelerate down- 
ward. This fact was supported by simulation runs which showed that it takes 
much more time to decrease y than to increase it. In effect, if y must be de- 
creased, then it varies more slowly than 3. The second point we wish to make 
is that time scale separation is also determined by sensitivity of the long term 
cost-to-go. Thus, it is possible that 3 should be considered slower than its 
ranking in equations (1.21) - (1.26) by virtue of its effect on the long range 
cost-to-go. Unfortunately, we have no method for analyzing the sensitivity of 
the cost-to-go for such a complex optimization problem as the one at hand. 
Therefore, we considered both the conventional orderings (x,y), (E,3)» h, y and 
(x,y), E,-3, h, y as well as the ordering (x,y), E, h, y, 3 indicated by the 
scaling and our discussion given above. For the initial conditions studied we 
found that the latter ordering yielded faster times than the conventional ones 
(mainly due to the asymmetric behavior of y dependent on whether it is increasing 
or decreasing). Note that if 3 is considered to be the fastest variable, then 
the SPT algorithm is considerably simplified. One essentially calculates a and 
u for vertical plane flight and adjusts a continuously to head to the target 
(or the predicted position of the target). The other orderings are considered 
in Section 5.4 (E faster than 3, 3 faster than h or y) and Section 5.5 (E,3 
together faster than h or y). 

In the following sections we will first work through the analysis required 
for the solution of each of the control subproblems, and then we will summarize 
the resulting algorithms and note the computation required for their on-board 
execution. As we noted in Section 4.2, the SPT approximation involves solving 
a hierarchy of control problems which must be solved in sequence (see Figure 
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4.2.2). Thus, the total time required to calculate one update of the aircraft 
controls by means of the SPT method is the sum of the computation times required 
for each of the algorithms which solve the individual subproblems. In addition 
to noting the computation time required for the individual subproblem algorithms, 
we will also note which functions must be stored and which quantities could be 
obtained directly from aircraft sensors. 

The organization of this chapter roughly follows the separation of time 
scales we have assumed for the aircraft model. Thus, Section 5.2 discusses the 
x,y reduced order problem. Section 5.3 discusses the E boundary layer problem, 
and Section 5.4 discusses the 3 boundary layer problem. Due to the difficulty 
of computing the 3 boundary layer solution we discuss the use of a suboptimal 
solution in Section 5.5. Then in Section 5.6 we discuss the h boundary layer 
problem and finally, in Section 5.7 we discuss the y boundary layer problem which 
operates on the fastest time scale. 

In each section we have discussed the computational requirements as well 
as the analysis required for the solution. In addition, in Sections 5.2 and 
5.3 we have discussed the possibility of linearization to solve the faster 
boundary layer problems corresponding to the subproblems of those respective 
sections. Thus, in Section 5.2 we consider the possibility of linearizing around 
the x,y reduced order solution as a nominal value to obtain linear cotrols for 
the E,3>h,y boundary layer problem. Likewise, in Section 5.3 we consider lin- 
earizing around our suboptimal solution for E to obtain linear controls for the 
corresponding faster boundary layer 3» h and y. 

Figure 5.1.1 shows the complete control logic for. long range interception. 

It involves an iterative calculation of the intercept point and the intercept 
time, t.p. In addition, the calculation of heading on the cruise arc and controls 
for the three parts of the trajectory (before cruise, cruise and after cruise) 
is also performed. 

5. 2 The Reduced Order Problem in x,y: The Cruise Solution 

5.2.1 x,y Reduced Order Control Problem and Feedback Solution 

By considering x and y on the same time scale and all the other states in 
(1.3) - (1.8) on a faster time scale we obtain the reduced order dynamical system 
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Figure 5.1. 1 

Overall Control Logic for Long Range Interception 
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(2.1) 

dx 

^ = V cos3cosy 

(2.2) 

^ = V sin3cosy 

(2.3) 

mg 

(2.4) 

0 _ 0- + uT^^^sina)sina 
^ mV cosy 

(2.5) 

0 = V siny 

(2.6) 

(L + uT„^sinc3t)cosa - 1 
A . njflx 

° iS 


together with the constraints (1.9) - (1.13). The Hamiltonian for this problem 
is 

(2.7) H = 1 + X V COS3COSY+X V sin3cosy 
A y 

and since the problem is to minimize time, the Hamiltonian is identically 0 along 
the optimal trajectory. The pseudocontrols E,h, 3, y and the controls a,a,u are 
chosen to maximize the velocity V subject to the equality and inequality con- 
straints and satisfy the x,y terminal conditions. Thus, one easily finds that 
the optimal pseudocontrols for the reduced order problem are given by: 


(2.8) 

E* = 


(2.9) 

3* = 

-1 

X -X 
Xf x^ 

(2.10) 

h* = 

"c 

(2.11) 

Y* = 

0 

(2.12) 

a* = 

“c 

(2.13) 

a* - 

0 

(2.14) 

u* = 

“r 
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where the subscript "c" denotes constant cruise values. Let denote the cor- 
responding maximum cruise velocity determined from and h^ by (1.18). For 
our particular aircraft, which is taken from Parsons (1972), these constants 
are given by 

(2.15) = 2.9949 X m 

(2.16) = 1.2192x10^ m 

(2.17) = 1.509° 

(2.18) = .8928 

(2.19) = 590.2 ms"^ 


In order to define the cost criterion for the E boundary layer we need to 

obtain the adjoints X and X for this reduced order problem. This is easy to 

A y 

do and we find that they are given by the expressions 


( 2 . 20 ) 

(2.21) Xy 


COSg^ 

sin$^ 


where we have introduced the notation 6^ for the angle given by (2.9) or by 

X “ X • 

(2.22) cosB. = . 5-jr- 

((Xf-x^)2+(y^-y.)2)’5 


yf-yi 

(2,23) sin3^ = « n-p 

"" {(x^:-x^)2+(y^.y.)2)^a 


Note that (2.22) and (2.23) uniquely define 3^ as the heading angle between the 
initial point (x. ,y..) and the final point (x.p,y^), measured from the positive 
X direction to the line between these points. Note that we can rewrite the 
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control law for 3* in (2.9) as 

(2.24) 3* = 3^ 

For convenience we will use the simpler expression (2.24) in the future rather 
than the expression (2.9). However, it should be kept in mind that 3^ is a 
function of the current x and y states as expressed in (2.9) or in (2.22), (2.23) 

The cost criterion for the next level of the SPT hierarchy of control sub- 

problems is found by substituting the expressions for the adjoints^ and A^ 
from (2.20), (2.21) into (2.7) to obtain 

(2.25) L = 1 - X cosycos(3-3-) 

V ^ 

5.2.2 Computational Requirement for the Solution 

Clearly, the computations required in this reduced order problem are trivial. 
All the controls and pseudocontrols are constants which are precomputed and 
stored, except for the pseudocontrol 3*=3^,, which is easily computed using one 
of the inverse trigonometric functions from (2.9), (2.22) or (2.23). Figure 
5.2.1 summarizes the computational requirements at this level. 

Storage 

(constant) 
h^ (constant) 

(constant) 
u^ (constant) 

Figure 5.2.1 

Computational Requirements of x,y Reduced Order Problem 

Note that these are the minimum on-board calculations that must be made and the 
minimum storage required for the x,y reduced order problem. At faster levels 
we will assume that E , h , a , u have already been stored and that 3« has 

^ w ^ ^ 

already been calculated, and we will not count this storage or calculation in 
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Calculation 

(by means of 
equation (2.9) 
or (2.22) or 
(2.23)) 


the computational requirement for the faster level solution. In this way, the 
total on-board computational requirement for the problem will be obtained by 
adding together the requirements at each level without the risk of counting 
storage or calculations more than once. 

5.2.3 L inearization Arou nd the x,y Reduced Order Solution > ' '■ 

In Subsectibn 5.3.4 we discussed the possibility of linearizing the boundary 
layer problem around the reduced order solution. In this case we would obtain 
a linear, quadratic criterion, time invariant, infinite time problem in the per-'' 
turbations of E, B, h, y> u, a and a from their optimal values computed in the ' 
x,y reduced order solution. The linear control problem has a simple control law 
of the form 



where G is a constant 3x4 gain matrix computed off-line. 

The equation (2.26) gives us the linear regulator control to maintain the , 
cruise conditions of the aircraft. Thus, on the cruise portion of the trajectory 
the full control law for the aircraft is particularly simple. Figure 5.2.2 
summarizes the on-board computational requirements in this case. 


Storage 


Calculation 

G 


u, a, a 

(a 3 x4 


(by means of 

array of 
constants) 


equation (2.26)) 

) 1 

Figure 5.2.2 



On-Board Computational Requirements During Cruise 
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5. 3 The Initial Boundary Layer Problem in E 


5.3.1 E-Boundary Layer Problem and Feedback Solution 

By considering E on a faster time scale than x or y and on a slower time 
scale than B, h or Yj we obtain the following boundary layer dynamical system 
for E: 


(3.1) 

(3.2) 


dt mg 

(L + uT ^ sina)sina 
g _ j max ^ 

mV cosy 


(3.3) 

(3.4) 


0 = V siny 


0 = 


(L + uTj^^^sina)cosa - mg cosy 
mV 


together with the constraints (1.9)- (1.13). The cost criterion for this probl 
is found from (2.25) to be 


(3.5) J = / [1 - V ^osycos(B-Bp)] dr 

The initial condition is given by E = but the final condition is given by 
E = E^ and not by E = E^. In particular, we see that the resul ting^ SPT control 
law will fly the plane to the cruise condition but will not yield a control 
which will give some final energy other than the cruise energy. We have dis- 
cussed this problem more generally in Subsection 4.4.2. We discuss a possible 
solution for this particular aircraft problem in Section 5.5 of this chapter. 

Since the control problem for this boundary layer has only one state dimen 
sion, we may use the technique of Subsection 4.3.2 to solve for tf?e feedback 
control law by minimizing or maximizing the expression 


(3.6) 


mg(l - Y cosycos(B-B^)) 

V(uT ^^cosa - D) 

^ max ’ 


with respect to the pseudocontrols B,h,y and the controls u,a,ci such that the 


Inequality constraints of (1.9)- (1.13) and the equality constraints of (3.2)- 

(3.4) are maintained. In addition, if we are minimizing (3.6), then we must 

restrict ourselves to pseudocontrol and control values such that ^ given by 

(3.1) is positive, and if we are maximizing (3.6), then we must restrict our- 
dE 

selves such that ^ is negative. 

From the equality constraints it is not hard to deduce that 

(3.7) Y = 0 

(3.8) 0=0 

Thus, the pseudocontrol y and the control a are unchanged from their values given 
by the x,y reduced order solution. Moreover, it is not hard to see that the 
optimization problem gives 

(3.9) B - 6^, 

so that the pseudocontrol B also has its reduced order value. Substituting (3.7), 
(3.8) and (3.9) into the original E-boundary layer problem, we obtain the fol- 
lowing optimization problem. We must maximize or minimize the function 

mg(l "Y ) 

(3.10) -£ — - 

subject to the equality constraint 

(3.11) 0 = L + uTj^gj^sina - mg 

together with the inequality constraints (1.9), (1.10), (1.11), (1.12). Using 
the equality constraint (3.11) to solve for u and substituting this expression 
back into (3.10) we obtain the expression 

mg(l-Y ) 

(3.12) 

V([mg - L]cota - D) 

which we must maximize or minimize, according to the convention described above 
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after (3.6). By eliminating u we also obtain a new inequal ity constraint on V 
and a corresponding to the original constraint on u in (1.12). The new constraint 
is 


(3.13) 


0 ‘< 


mg •• L 
max 


ke have now reduced the optimization problem to a problem involving only 
two variables, V and a. Let us now consider V fixed and examine the dependence 
of (3.12) on a. Using (1.14)- (1.16) we can write the a dependence explicitly 
in (3.12) as 


(3.14) 


mg(l-^ ) 
c 


V([mg- L^alcota- Dq - L^na ) 


Similarly, the inequality constraint (3.13) becomes 


(3.15) 0< 


mg - L^a 


a 


max 


< 1 


It is not hard to see that for a in the range 0<a<7r the denominator of (3.14), 
namely 

(3.16) V([mg - L^a]cota - Dq - L^na^) 

is a strictly decreasing function of a. Thus, the expression (3.14) is an in- 
creasing function of a in this range. Note that this range easily includes all 
admissible a since cXg < 'tt (in our aircraft example 

Since (3.14) is monotonical ly increasing in a, in order to maximize (3.14) 
with respect to a we must choose the largest possible a consistent with the in- 
equality constraints (3.15) and (1.9). Similarly, if we wish to minimize (3.14), 
we must choose the smallest possible a consistent with the inequality constraints. 
Note that the function 
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mg - L a 

(3.17) u(a) = Y “ 


max 


sina 


which appears in the inequality constraint (3.15) is monotonically decreasing in 
CL for admissible a (in fact, for all a such that 0<ot<^). The mbnbtonic beha- 
vior of the function (3.17) allows us to describe the maximum and minimum allow- 
able values of a fairly easily. There are essentially three different cases to 
consider. 

Case 1 . This case corresponds to 
(3.18) L^«s + W’"“s ^ '"9 


In this case there is no admissible a which satisfies the inequality constraints 
(3.15) and ’(1.9). Physically, (3.18) represents the stall constraint for the 
plane and it defines the region in which the plane cannot maintain vertical equi- 
librium using both lift and thrust forces. 

Case 2 . This case corresponds to the two relations 
(3.19) < mg 

and 

<^•20) - ""S 

Physically, there is less lift force than the weight of the aircraft, but by 
using maximum thrust in addition to lift, the plane is able to maintain equili- 
brium. The maximum admissible a is given by a = a^. The minimum admissible a is 
given by solving the equation 

(3.21) " "^9 

Case 3 . This case corresponds to 

(3.22) > mg 

Physically, the aircraft has enough lift to maintain vertical equilibrium by 
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lift force alone. The maximum admissible a is given by 

(3.23) a = ^ 

a ' . 

and the minimum admissible a is given by solving (3.21) as above. 

The argument above concerning a allows us to express a as a function of V 
and E, and hence allows us to eliminate a from the problem to optimize (3.14).' 
Since the aerodynamic quantities L^, Dq and T^^^^ are given only as data, we must 
solve the V optimization numerically. We have done this off-line and stored the 
resulting optimal control laws V*(E). Note that there are two V* controls, one 
corresponding to the maximization and one corresponding to the minimization of 
(3.14). We have found that for the optimal V = V*(E), only Case 3 above occurs. 
Thus, the throttle control u is 0 or 1 depending on whether E is above or below 
the cruise energy E^. The angle-of-attack a is chosen to maintain vertical equi- 
librium, i.e,, to satisfy equation (3.23) if u = 0 and to satisfy (3.21) if u = l. 
In Figure 5.3.1 we have summarized the SPT feedback control law for the E-boun- 
dary layer. 

To obtain the cost criterion for the next level control problem, we need 
to obtain the optimal adjoint A| for the E boundary layer 'problem. This is 
easily accomplished by solving 

(3.24) 0 = (T*cosa*-D*) 

Note that (3.24) gives A| as a function of V*, which in turn is a function of E, 
namely 


V* 

mg(-y — 1) 

(3.25) A* = — ^ 

'' V*(T*cosa*- D*) 

We may either store A| explicitly as a function of E or we may compute A| from 
V*(E) using (3.25). However, the latter approach requires us to store the aero- 
dynamic functions T*, D*. 
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5.3.2 Computational Requirements for the Solution 

The important new quantities computed in the E level control problem are 
the pseudocontrol s V*(E) and h*(E), the controls a*(E) and u*{E), and the adjoint 
X|(E). We have stored V* (for both E < E^ and E>E^) as a function of E. From 
V*(E) we can easily compute h*(E) as 

(3.26) h*(E) = E - 

The control function u*(E) is a trivial switch function as can be seen from 
Figure 5.3.1. For E>E , we have u*(E) = 0 and for E < we have u*(E) = 1. 

The control a*(E) is obtained from solving (3.21) if E<E^ and from solving 
(3.23) if E>E^. Equation (3.21) gives a only implicitly so that we must approxi- 
mate the exact solution. For L^» mg, the linear approximation sina-a gives an 
explicit formula, 

(3.27) a*(E) = j '^9 

max a 

for a*(E). Using Newton-Raphson' s method with (3.27) as an initial guess, we 
can obtain a better approximation of a*(E) with only a few iterations. 

On-board computation of a*(E) from either (3.21) or (3.23) will require 
storing the aerodynamic functions T^^^(E,h) and L^(E,h). The function has the 
form 

(3.28) L = ^.C. (M)p(h)V^S 

a 

where S is a constant, M is Mach number related to velocity by the equation 

(3.29) M = ^ 

where c(h) denotes the speed of sound at altitude h. Thus, to store L^(E,h) it 
suffices to store functions of only one variable and then compute from these 
using (3.28). 

The function T on the other hand, is Only given as data and is not 
given in terms of functions of a single variable. Thus, we must store T^^^ in 
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an efficient form. At present T has been approximated using bicubic splines, 

md A 

but the approximation is costly in terms of storage and this method of storing 
^max improved. 

Figure 5.3.2 summarizes the computational requirements at the E-boundary 
layer level . 


Storage 


Computation 


V*(E) (function of 
one variable) 
(piecewise linear 
approximation over 
five intervals in E) 

L (E,h) (function of 

two variables)t 

(function of 
two variables) 


a* 

(from (3.21) 
or (3.23)) 


h* 

(from (3.26)) 


Figure 5.3. 2. -Computation and Storage Requirements for E-Layer Solution. 


5.3.3 Linearization Around the E-Boundary Layer Sol u tion 

As in Subsections 4.3.4 and 5.2.3, we can also linearize the boundary layer 
control problem for B, h and y around the nominal values for these variables 
obtained in the E-boundary layer problem solution. As before, we obtain a linear, 
quadratic cost criterion, time invariant, infinite time problem in the pertur- 
bations of the variables B, h, y» ot from their optimal values computed in the 
E-boundary layer problem. Note that we assume that u retains the same value, 
either 0 or 1, in the higher level control problems that it has in the E level 
problem. The linear control problem has a simple control law of the form 

(3-30) (0 = (a*(E) 

where G is a 2x3 gain matrix which depends only on E. In this case the gain 
matrix G is a function of the energy E and must be computed off-line and stored 

tAlternati vely, store the functions C, ,p,c of one variable and compute L from 
(3.28). a “ 
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efficiently for on-board use. Note that storage of G amounts to the storage of 
six functions of one variable. Figure 5.3.3 summarizes the computational re- 
quirements for this linearization. 


Storage 

G 

(a 2 X 3 array 
of functions 
of one variable) 


Calculation 
a, a 

(by means of 
equation (3.30)) 


Figure 5. 3. 3. -Computation and Storage Requirements for Linearized 
E- Layer Solution. 


5.4 The Initial Boundary Layer Problem in g 

5.4.1 g-Boundary Layer Problem and Feedback Solution 

By considering g on a faster time scale than x, y or E and on a slower time 
scale than h or we obtain the following boundary layer dynamical system for g: 


(4.1) 

(4.2) 

(4.3) 


dB _ (l- + uWina)sino 
dt mV cosy 

0 = V siny 


0 = 


(L + uT^^^sina)cosa - mg cosy 
mV 


together with the constraints (1.9) -(1.13). The cost criterion for this problem 
is found to be 


^00 

(4.4) " Y ^°^YCOs(g-g^) + a* ^ (uTj^g^cosa - D)] dx 

where A| is the optimal adjoint computed in (3.25). The initial condition is 
given by g=g^. but the final condition is 3=3^.. Again, as for the E-boundary 
layer problem, we see that the SPT control law will fly the plane to the cruise 
heading but will not yield a control which will give the final heading angle. 

The optimal control law is obtained by minimizing or maximizing the expres- 
sion 
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[1 - Y cosycos(6-B^.) + 

(4.5) ^ — — 

(L + uTj^g^j^sina)sina 

with respect to the pseudocontrols h, y and the controls u, a, a such that the 
inequality constraints of (1.9) - (1.13) and the equality constraints (4.2) and 

(4.3) are satisfied. If B < 3_ then we minimize (4.5) and we must restrict our- 

^ ds 

selves to pseudocontrol and control values which give a positive value for ^ 

in (4.1). Likewise, if 3 > then we maximize (4.5) and we must restrict our- 

C JO 

selves to pseudocontrols and controls such that ^ remains negative. 

From the equality constraint (4.2) it is not hard to see that 

(4.6) Y = 0 

At this stage we will make the additional assumption that u takes on the same 
value (namely 1 or 0) that it does in the E-boundary layer solution. In general, 
we assume that if a control variable takes only a maximum or minimum value at a 
slower level of the SPT hierarchy, then it maintains that value of the control 
at the faster levels of the hierarchy. Thus, instead of writing uT^,„ we will 

maX 

write T and assume that u is already determined and fixed by the E-boundary 
layer solution.* 

With these assumptions the expression (4.3) becomes 

(4.7) 0 - (L + T sina)cosa-mg 

Using (4.7) to simplify the denominator of (4.5), we obtain 


(4.8) 


[1 - Y cos(3-3^,) + 

g tana 


Let us fix V and try to determine the optimum a and o for a given fixed V. To 
do this we consider a a function of a determined by the expression (4.7). Taking 

*A comparison with, the exact energy state solution using both E and 3 as states 
on the same time scale shows that this approximation is valid in all cases 
except those where level flight cannot be maintained during a turn with zero 
thrust (see Figure A5.1.9). This special case is handled quite easily in 
practice by switching to full thrust if level flight cannot be maintained other- 
wise. 
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the derivative of (4.8) with respect to a, we obtain 
-A*V 

(4.9) -jj|- (T sina+2nL^a) I cota 

- : - (1 - ilfCT cosa - - nL^a^]) | sin’^a |f 

The derivative is determined from (4.7) to be 


(4.10) 



L + T cosa 

— cota 

L a + T sina 
a 


If we assume that a is small and that we can neglect terms higher than linear 
terms in a, then we can determine the optimum a by setting the expression (4.9) 
equal to 0, In this way, we obtain 


(4,11) 


2 

tan a 


1 - ^^cos(3-3^) + (T - Dq) 


«X|V 


(T + 2nL ) 


a' 


(T + L y 
a 


mg 


In (4.11) one takes the positive or negative root depending on whether 3 < 3^ or 
3>3^* The corresponding value for a is determined from (4.7) using the value 
for a from (4.11). Note that we must also test the possibility that a or a lie 
on a constraint. One does this simply by setting solving (4.7) for 

the corresponding a or by setting a = a^ and solving (4.7) for a.' We then test 
these values against the value of a determined by (4.11) and keep the one which 
gives the smallest (if 3<3p) or greatest (if 3>3_) value of (4.8). 

So far we have only determined the optimum a and a values in terms of a 
fixed V. To obtain the optimum V, and thus solve the problem, we must resort to 
numerical optimization off-line.. This optimization is made difficult by the 
dependence of (4.8) on E and 3”3^, Thus, the computed optimal v-elocity pseudo- 
control^ V* for the 3-boundary layer is a function of the following form 

(4.12) V* = V*(E,3-3^.) 
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Note that V*(E.3-6^) = V*(E,-(3-6^) ) . 

The control law for the 3-boundary layer is suninarlzed in Figure 5.4.1. 

5.4.2 Computational Requirement for the Solution 

The major computational requirement at this level is the storage of 
V*(E,3-3j,)« Having stored V*, we can compute o* from (4.11) or from the con- 
straint values and then compute a* from the expression 

(4.13) a* ^ a Jp* 

max 

We will also have to compute h* corresponding to V*(E,3-6j,) from (3.26), and we 
need X|(E) to obtain a* from (4.11). If is not stored, then it must be com- 
puted in the E-boundary layer by means of equation (3.25). Note that in (3,25) 
the expression V* is V*(E), the optimum V for the E boundary layer problem, and 
not V*(E,3-3^). The computational requirements at the 3-boundary layer level 
are summarized in Figure 5.4.2. 


Storage 


Calculation 


V*(E,3-3^,) (function of 

two variables) 

Dq (function of two 
variables)t 

n (function of one 
variable) 


(from (3.25)) 
a* 

(from (4.11)) 
a* 

(from (4.13)) 


Figure 5.4.2 

On-Line Computational Requirements for the 3-Boundary Layer 


In the next section we take an alternative, simpler approach to the 3- 
boundary layer solution discussed here. We use the work of Parsons (1972) to 
obtain an exact solution to the reduced order problem in x,y,E,3. Then we ap- 
proximate this exact solution with a suboptimal solution that essentially divides 
the E-3 plane into regions where zero or maximum thrust is used and all the 

tAl ternatively, store functions C^ ,p,c of one variable and compute Dq. 
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Figure 5.4.1 Boundary Layer Algorithm, 
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turning takes place either on the maximum turn rate (MTR) locus or on the maxi- 
mum velocity constraint boundary. This approach involves more storage, but less 
on-line computation. 

Linearization around the 3-boundary layer solution for h and y can be per- 
formed in the same way as the linearization around the E-boundary layer (see 
Subsection 5.3.3). 


5.5 E,3 Boundary Layer Problem 
5*5.1 Exact E,3 Solution 

Due to the difficulties of computing the 3 boundary layer feedback control 
described in Section 5.4, we now suggest an alternative approach. The work of 
Parsons (1972) provides us with an exact solution to the x,y,E,3 problem in 
which x, y, E and 3 are treated on the same time scale. Using this exact solu- 
tion as a basis, we develop a suboptimal feedback law which can be computed on- 
line and which is not far from optimal. 

Assuming that h and y vary faster than x, y, E, 3 (which vary on the same 
time scale), we obtain from (1.3) -(1.8) the dynamic equations 


(5.1) 

(5.2) 

(5.3) 


d Y 

^ = V cos3cosy 
^ = V sin3cosy 


dE _ 
dt 


mg 


(L + uT^vSi'>“)si 


^ '"max 

dt 


na 


(5 4) M =: 

' * ' dt my cosy 

with the equality constraints 
(5.5) 0 = V siny 

(L + uT„,„sina)cosa - mg cosy 

niaX 


(5.6) 


0 = 


"inV* 


from (5.5) we have y =0 and from (5. 6) we| obtain 
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(5.7) mg seca = L + uT|^^j^sina 


Substituting this expression in (5.4), we have 


(5.8) 


d3 _ q tana 
'St V 


Using the fact that a is small and L is much greater than we can find 

a max 

an approximate expression for ot from (5.7), namely 


(5.9) a = 


Substituting this into (5.3) and approximating uT^^^cosa by uT we obtain 

maX iTluX 


(5.10) - 


mg 


where D|^ is given by 

/ 

(5.11) 


Thus, the 

(5.12) 

(5.13) 

(5.14) 

(5.15) 


system (5.1) -(5.4) becomes 




^ = V sine 


dE . ''(uT^ax-‘’0-‘’L=®'= 


dt mg 

d3 _ g tana 
dt ■ V 


There are inequality constraints on a, directly through (1.11) and indirectly 
through (1.9) and (5.9). 

This optimal control problem (minimum time) must be solved numerically in 
terms of a TPBVP derived from the first order necessary conditions. In Appen- 
dix 5.1 we set up the TPBVP and discuss the optimization procedure. In the next 
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section we describe the suboptimal approximation we have chosen for this solu- 
tion. 


5.5.2 Suboptim a l Real-Time Approximation 

This subsection deals with the problem of computing in real time a subop- 
timal minimum-time trajectory to intercept a moving object, by approximating 
the work of Appendix 5.1. These sections have shown that optimization on-line 
is too time-consuming even for the energy-state approximation and therefore not 
feasible in real time. A sacrifice in accuracy is called for so that a near- 
optimal solution may be generated on-line merely by looking up stored trajectories 
and fitting together portions of them to match the particular boundary conditions 
imposed upon the aircraft states. The program RLTIME has been written to produce 
such a near-optimal trajectory. Boundary layers can then be added to this 
solution in order to smooth out the sudden changes in the states of the aircraft 
introduced by the energy-state approximation. 

Stored Trajectories . Four trajectories are stored: 

(1) Max -turn locus with 9022 <E< 29870 m. This locus has been defined in 
Appendix 5. l--basically it requires a =12° and thrust is maintained 

either at T = T^,„ or at T = 0. The upper limit on E was chosen arbitrarily while 
the lower limit is approximately the Lufbery Circle Point value (Parsons (1972)). 
This is the classical aerial combat situation and is a steady-state turning 
rate condition. If the airplane were to fly along this locus its energy would 
decrease while following a spiral -shaped trajectory in three-dimension space. 

The following equations were integrated to obtain the locus with T = T^g^^:> 


(5.16) 

(5.17) 

(5.18) 


ds - 

V 

V„t(T„^„ - D„ - D, sec^CT ) 
max 0 L m' 

M cosB 

(T - D„- D, sec^a ) 

' max 0 L m 


dt _ 
^ ' 


dx 

dE 
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(5.19) ^ = 


W sin3 


(^max - Dq - %'> 


Equations (5.16) - (5.19) have been derived from the equations of motion of the 
airplane model (5.12) - (5.15). Note that for 9022 < E< 29870 m, ^ DRAG and 
so for ^ increases as energy decreases. Vj^j is the velocity necessary 

to remain on the locus and is given by the equation: 

(5-20) 

\ a 



which merely states that LIFT=W seca at maximum bank angle and angle-of-attack, 
i.e. , setting = a^. 

In order that the x and y values obtained by integration of equations (5.18) 
and (5.19) may be translated into the inertial coordinate system* consider the 
following: 


f(E) from (5.16), hence 


(5.21) 


3(E^) - 3(Eq) = 31^(5^) ” 



f(E) dE 


where the subscript MT refers to the max-turn locus coordinate system and the 
unsubscripted variables refer to the inertial coordinate system. (See Figure 
5.5.1.) Equation (5.18) may now be rewritten as: 


(5.22) 

(5.23) 


dx 

dE 


W cos[3^y(E) + (3 (£q) “%(Eo))1 


(T-D) 


“ cos{3(Eq) - 


-W cos3j^(E) 


- sin(3(EQ) - Bj^-j-(Eq)) 


T- D 

-W sin3j^-p(E) 
T^Td 


Hence, 
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(5.24) x(E) - x{Eq) = cos(6(Eq)-B„^(Eq))[x„t(E)-x„^(Eq)] 

+ sin(B(EQ) - B|^j(EQ))[y^j(E) -yj,|j{Eg)] 

and similarly: 

(5.25) y(E) - y(EQ) = cos($(Eq) - 3^t(Eq) )[ y^y(E) - y^-^(EQ)] 

+ sin(B(EQ) - Bj^-|-(Eq) )[X j^y(E) - X|^j(Eq)] 

(2) Max-turn locus with 8992<E< 1524 m. Here the lower limit is arbitrary 
while the upper limit is again approximately the Lufbery Circle Point value. 

To obtain the turning in the same sense as in (1) with a + ve, it is necessary 
to move up the locus with increasing energy since now E^Lisitions 

(5. 16) - (5. 19) were integrated and again equations (5,24) and (5,25) give the 
horizontal plane distances in the inertial frame, 

(3) The AB=0> minimum-time energy path to the cruise arc . This is a 
trajectory belonging to the family of trajectories presented as initial turns 

to the cruise arc in Appendix 5.1. This particular one allows no turning and is 
really just the vertical plane solution of the problem of reaching the cruise 
arc in minimum time. 

(4) The AB=0 chatter path off the cruise arc . Appendix 5.1 presents the 
theory supporting this trajectory, which requires zero thrust and is a maximum 
deceleration arc. 

Real-Time Approximations . As stated earlier, a real-time solution requires 
minimum computation on-line and should involve mainly looking up of stored tra- 
jectories and fitting together portions of them so as to satisfy boundary con- 
ditions. One possible approximation for real-time determination of trajectories 
is to turn only on the max-turn locus and at other times to follow either the 
delta-beta =0 path to the cruise arc or the chatter path off the cruise arc. 
Hence, for initial turns to the cruise arc, the approximate solution requires 
a zoom to the max-turn locus from the initial conditions followed by turning 
on the locus and then a zoom to the delta-beta = 0 path to the cruise arc. 

As presented in Appendix 5.1, Figure A5.1.9 shows the optimal trajectories 
in the beta-energy plane for final turns from the cruise arc. It is possible 
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to make approximations in varying degrees to produce near-optimal trajectories 
to replace the optimal ones, in real time. Using the same approximation as for 
the initial turns, all the final turning is also on the locus. Thus the real- 
time approximation requires bank angle chatter on getting off the cruise arc 
followed by a zoom to the max-turn locus. After turning on the locus a zoom at 
constant energy to the final conditions completes the trajectory. This approxi- 
mation is indicated by the dotted lines in Figure 5.5.2 where the thrust-switch 
locus is also ignored: chatter is over. 

Numerical results have been obtained for the approximation described above. 
Probably, a better approximation would be as in Figure 5.5.3 where a piecewise 
linear fit is made for the thrust-switch locus and the zoom-climb locus and the 
trajectory allows turning before reaching the max-turn locus by about 60°. The 

rate of turning may be approximated as constant by the almost straight line tra- 

dE 

jectories, or al ternatively a crude ^ integration may be performed on-line. 
Switching the thrust would require storage of the max-turn locus with thrust=0. 

Formulation of the Problem . We assume for simplicity that the target is 
moving in a straight line at constant velocity. For interception of the target, 
we require pursuer trajectory time (tj^) = target time to interception point (t^). 
In addition t^ = t^ must be the minimum possible time. Hence an iterative loop 
is necessary to match up t^ and t^ for minimum-time interception. With a value 
for t^, the position of the target in three-dimensional space may be obtained 
and the problem may be restated as: given initial and final energies, altitude, 

velocity, and horizontal plane positions, a minimum-time trajectory is required 
so that the above boundary conditions are met. 

Figure 5.5.4 shows how the min-time path is split up among the stored tra- 
jectories and the cruise arc. The heading angle, 3> is measured clockwise from 
the x-axis. Two constraints have to be observed. Restricting ourselves for 
simplicity to the case of clockwise turns, (1) the relative angle between the 
initial and final points (3^) is assumed greater than 6 q, the initial heading 
angle; (2) the total turning angle is assumed less than 180° 

In general terms the algorithm is: 

(1) Estimate t^, the interception time 

(2) compute final conditions for interceptor aircraft 

(3) from initial altitude, execute a velocity (or altitude) zoom (constant 
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VHO 


Suboptimal approximation (linear fit for thrust switch and zoom climb loci) 



Figure 5.5.3 Change in Heading for Variable-Altitude 

Final Turns from Cruise Arc 


.Pk 



Figure 5. 5. 4. -Geometry for Determining Initial Heading Angle and 
Interception Time t^. 


energy) to the max- turn locus 

(4) determine cruise heading 3^ (see equation (5.28) below) 

(5) turn along the max-turn locus until heading = 8^ 

(6) zoom to the A3=0, min-time E-path to cruise arc 

(7) follow A3=0, min-time E-path to cruise arc 

(8) stay on the cruise arc for distance 

(9) follow A3= 0 chatter path off the cruise arc for distance R 2 

(10) zoom dive or climb to max-turn locus 

(11) follow max-turn locus until heading angle = 3^ 

(12) zoom dive or climb to required final altitude and velocity 

(13) compute trajectory time and update t^ using Newton-type step 

(14) go to (2) if accuracy on t^ not met. 

Figure 5.5.8 introduces relationships between the various distances involved 

(5.26) >'1 + ^2 ^^l"^ R 2 + R^)sin3^ ^ ^f “ ^0 

(5.27) ^2 ^f " ^0 

Which give the cruise heading as: 



and since Xj, y|» ^^ 2 * ^2 functions of 3^> equation (5,28) has to be solved 
iteratively, for example by using a Newton method. R^ is then obtained from 
equation (5.26) or (5.27). 

5.5.3 Computational Requirements 

The two important quantities being calculated in this algorithm are 3^ and 
t^, both by iterative methods involving just a few iterations each. We consider 
here computational requirements for calculation of each of the above two quan- 
tities which are s( fficient for the second (more accurate) real-time approxima- 
tion mentioned in the last section. Most of the trajectories to be stored may 
be approximated by piecewise linear fits to keep storage requirements dov/n. 
Figure 5.5.5 summarizes the computational requirements. 

In the example case presented and considered in Chapter 6, 3^. = 38.3° while 
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storage 


T = 0 1 Ax(3,E) 

T = j^y(6.E) 

on max-turn locus 


Computation 


using 

equation (5.28) 
in iterative loop 


thrust switch locus (B,E) 


T = 0 

. At(E), h(E) 

^ ^max 

on max-turn locus 


thrust switch locus (BjE) 
zoom climb locus (B^E) 

3-E lines on max-vel 
constraint (see Figure 
5.5.2) 


iteratively using 
t^ = sum of times on 

separate portions of 
trajectory; t^ time 

R 

on cruise ~ 5 g § ~2 


h(E), At(E), ax(E) on 
A 3 = 0 path to cruise 


h(E), At(E), ax(E) on 
AB = 0 chatter 1 ine 
(max-vel constraint) 


Figure 5. 5. 5. -Storage and Computation Requirements for E,B Layer. 



the relative angle between the. initial and final points was 38.65°. Hence, the 
3^ calculation may be omitted to save computation time with little loss in 
accuracy. 


5. 6 Initial Bo undary Layer Problem in h 

5.6.1 h- Bound ary Layer P robl e m a nd Feedback Solution 

As we discussed previously in Sections 5.4 and 5.5, we have separated the 
vertical plane part of the trajectory from the horizontal turn part. Thus, in . 
considering the h-boundary layer we need only look at the vertical plane case 
with a=0 and 3=3^. Moreover, we assume that any controls which take their 
value on a constraint in slower levels retain the constraint value at faster 
levels.* Consequently, we assume that T is determined (u = 0 if E > and u = l 
if E<E^). With these assumptions the dynamical system for the h boundary layer 
becomes 

(6.1) ^=VsiriY 

with the equality constraint 

(6.2) 0= (L + T sina)-mgcoSY 
and the cost criterion 

^00 

(6.3) J =J^ [1 - Y C0SY+ A* ~ (T cosa- D)] dx 

where A| is the optimal adjoint computed in (3.25). The initial condition is 
h == h^. and the final condition is h=h*(E), where h* is determined from (3.26). 

The optimal control law is obtained by minimizing or maximizing the expres- 
sion 

*More general solutions without these assumptions are given in the appendices. 

In particular, it is shown that, for turning maneuvers with a?^0, the results 
of Sections 5.6 and 5.7 are extended quite easily. 
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(6.4) 


[1 - Y COSY+ A* ^ (T cosa - Dq - nL^a^)] 

V siny 

with respect to the pseudocontrol y and the control a such that the equality 
constraint (6.2) is satisfied. If h<h*(E), then we minimize (6.4) and we must 
restrict ourselves to y>0. If h>h*(E), then we maximize (6.4) and we must 
restrict ourselves to y<0. Note that the expression (6.4) is an odd function 
of Y and that the expression (6.2) is an even function of y* It follows that 
the minimum solution y>0 is the negative of the maximum solution y<0 and that 
a is the same in both cases. Thus, we restrict our attention to the case when 
Y > 0. 

To determine the solution y>0 that minimizes the expression (6.4), define 
the Lagrangian function ^ as 

[1 - T7 cosy + Xt ~ (T cosa - - qL a^) ] , r, . t c-ir, mr, 1 

V E mg 0 ' a A[L a+T sina-mg cosyj 

(6.5) = ? — — + —^ 

V siny V siny 

where A is an undetermined Lagrange multiplier. Setting we obtain 

the relations 

(6.6) 0= (y +Amg^ “ cosa - Dq - a^nL^) + A(L^a + T si na)) cosy 

(6.7) ^ "^E sina+ 2qL^a) + A(L^ + T cosa) 

Assuming that a is small enough to neglect in (6.6), (6.7),* we find that Amg in 
(6.6) is negligible and we obtain the following expression for cosy 

V 

V, 

(6.8) ^ = cosy 

+ (T- Dn) 

E mg 0 

We solve (6.8) for cosy and then obtain a from (6.2). Note that for values of h 
*This approximation is quite good for needs refinement for T=0 case. 
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close to h*, the expression (6.8) may exceed 1 slightly (by .02) due to the 
absence of the nonzero a in the denominator. In this case we can either set 
Y*=0, since it will be very close to 0, or we can use a method of successive 
approximation to obtain a better approximation--for example, solve (6.2) with 
Y=0 and use this approximation of a in (6.6), (6.7) to obtain a better approxi- 
mation of Y- In Appendix 5.2 we have carried out a more careful analysis of 
the solution of equations (6.2), (6.6), (6.7) using perturbation analysis. Fig- 
ure 5.6.1 summarizes the control law computations for the h-boundary layer. 

The adjoint for the h-boundary layer is obtained from solving 

(6.9) 0=1--^ cosy+A|~(T cosa - Dq - nL^a^) + Aj!jV siny 

for A*, where y and a take their optimum values for the h-boundary layer problem. 
5.6.2 Computational Requirements 

The computational requirements at this stage are rather minimal. All storage 
requirements have been accounted for on slower levels. Likewise, the adjoint 
A| has already been calculated. Hence, the only computational requirement is to 
calculate y from (6.8) and a from (6.2). Figure 5.6.2 summarizes these minor 
requirements. 


5.7 Initial Boundary Layer Problem in y 

5.7.1 Y-Boundary Layer Problem and Feedback Solution 

We make the same assumptions for this problem that we did for the h-boundary 
layer. In particular, we assume that T is fixed, and a=0 and 3=3^.. The dynamic 
system for the Y'boundary layer becomes 

(7.1) sina-mg cosy 
with the cost criterion 

-CO 

(7.2) J ~J^ ” Y cosy+ a| ~ (T cosa - D) + Ajj^V siny] dx 
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Read in h 



Y* == Yn F" 


h < h* ( E ) h > h* ( E’) 

-J <J ? h*(E) > 


Figure 5.6. 1 .-h-Layer SPT Algorithm. 


Storage 


Calculation 


None 

(accounted for 
at slower levels) 


( from (6.8)) 


(from (6.2)) 


Figure 5.6. 2. -Storage and Computation Requirements for h-Layer 








The optimal control law is obtained by minimizing or maximizing the ex- 
pression 

1 - Y COSY+ ^ (T cosa - Dg - + Xj^V siny 

(7.3) — — 

L^a + T sina-mg cosy 

with respect to the control a. If Y<y*(E,h), then we minimize (7.3) and we 
must restrict ourselves to a such that + T sina-mg cosy>0. If y>y*(E,h), 
then we must maximize (7.3) and we must restrict ourselves to a such that 
L^a+T sina-mg cosy<0. 

Taking the derivative of (7.3) gives us 

(7.4) {-(L^a+T sina-mg cosy)(X| [T sina+ 2nL^a]) 

-(L^+T cosa) ( 1 - Y cosy + A| ^ (T cosa - D) + AjJjV siny)}(L^a + T sina-mg cosy) ^ 

When either E > E^ and y < y* or when E < E^ and y>y*, the expression (7.4) is 

negative. Hence, in these cases the expression (7.3) is decreasing in a. When 
E<E^ and y>y*, we want to maximize (7.3) so we choose a=0. When E > E^ and 
Y<y*, we want to minimize (7.3) so we choose a = a^. 

To solve for a in the other cases, we set (7.4) equal to 0 and we neglect 

terms in a higher than linear.* Define aQ by 

(7 5) a = fng. cosj]^ 
a 

Then setting (7.4) equal to 0 gives the following quadratic equation, 

(7.6) a^ - aQa + k = 0 
where k is given by 


*In Appendix 5.3 we have considered a higher order approximation. 
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(7.7) 


k = 


1 - Y COSY + ing + Aj^V siny 




Note that k<0 for E<E and k>0 for E>E_. The solution a is given by 


(7.8) a = ^ + !s-\^Q^-4k 
when E<Ej^, and 

(7.9) a - ^ - 4k 


when E > E^. Figure 5.7.1 summarizes the feedback law in this case. 


5.7.2 Computational Requirement 

As for the h-boundary layer, the computational requirements at this stage 
are rather minimal. All storage requirements have been accounted for on slower 
levels. Likewise, the adjoint has already been calculated. The adjoint 
must be calculated from (6.9) at this level. Finally, a is calculated from 
(7.8) or (7.9). The computational requirements are summarized in Figure 5.7.2. 


Storage 

None (all 
accounted for 
on slower 
level s) 


Calculation 
(from (6.9)) 

aQ 

(from (7.5)) 
k 

(from (7.7)) 
a 

(from (7.8) 
or (7.9)) 


Figure 5.7. 2. -Storage and Computation Requirements for y-Layer. 


A pictorial representation of the control problem is given in Figure 5.7.3. 
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Y-trajectories viewed on the 



Figure 5.7.3 

y-Boundary Layer Control 



The ytrajectories are shown on two time scales, t and T = t/e. The t time scale 

corresponds to events in the h-boundary layer and the change in y on this time 

scale appears instantaneous. Case I shows a trajectory for y*“ 0 aod y>y*> 

E<E^. From Figure 5.7.1, the optimal a=0 till y = Y* and a = [equation (7.8)] 

when Y<Y*. The asymptotic value of a is determined from the ot-value on the 

h-boundary layer, namely . Notice that a slight undershoot in y would 

‘"a 

occur and some form of anticipatory action may be used to reduce this undershoot 
and consequent chatter. Case II corresponds to an asymptotic value y = y*» 
y(t = 0)<y* and E’> E^. For this case. Figure 5.7.1 shows that the optimal 
value is till y(t)>y* and a= [equation (7.9)] for y<y*- In practice, it 
would be better to set a to its asymptotic value once |y-Y*| is less than a 
small threshold value. An alternative would be to switch to linearized control 
once |y-Y*l decreases below a prespecified value. A similar logic should also 
be used for other boundary layers. 
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APPENDIX 5.1 

DERIVATION OF THE EXACT E,3 SOLUTION 
A5.1.1 Introduction 

As we saw in Subsection 5.5.1 of Chapter 5, the equations of motion for the 
x,y,E,3 system (referred to as the energy-state approximation in Parsons (1972)) 
are 

. V(T - Do - D, sec^a) 

(1.1) E= 

(1.2) e = 9 tang 

(1.3) X = V COS3 

(1.4) y = V sin3 

E, 3» X and y are the four states and T, V and a are the three controls. The 
three states V, h and y of the original dynamic equations have been reduced to 
the single state E and V or h becomes a control variable. Essentially, the 
assumption introduced here is that V and h can be changed instantaneously at 
constant E by zoom climbs or zoom dives. For consistency with equations (1.3) 
and (1.4) these maneuvers also occur without changes in x or y. Section 5.6 of 
Chapter 5 presents the concept of boundary layer corrections to smooth out zoom 
climbs or dives as in Figure A5.1.1. 



Figure A5.1 .1 .-Concept of Boundary Layer Corrections. 
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The major assumptions implied in the equations of motion (1.1) to (1.4) 
are as follows: (a) change of mass is negligible, (b) cosine of angle-of-attack 

is approximately one, (c) the normat component of thrust ts negt-igihle compared 
to ti-ft. 

The boundary conditions that go with the formulation of (1.1) through (1.4) 
are for a turn to a point: 

(1.5) 

E(0) = Eg 

E(tf) = 

(1.6) 

8(0) = 

B(t.p) ~ 

(1.7) 

x(0) = Xq 

X(t.|r) = X^ 

(1.8) 

y(0) = yg 

y(tf) = y^ 

and for 

a turn to a line: 


(1.9) 

(1.10) 

8(tf) = tt/2 

X(tf) = Xf 

(rotate x-axis so that it is 
normal to objective line) 

(1.11) 

y(tf) = free 


The 

constraints on the three controls are: 

(1.12) 

0 = T < T < T„^^(h,M) 
nnn max ' 


(1.13) 

|al < a^(V,E,ag) 


(1.14) 

V < V (E) 
max'' ' 


(1.15) 

E 

b 

VI 



The minimum or idle thrust is .assumed to be zero, is the bank angle 
which yields the stall angle-of-attack = 12° at an existing energy and velocity 
from the equation: 

(1.16) = cos'^ f-j^] 
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/ 

J 

/ 

I 

/ 

j 

The V constraint combines an approximate maximum air speed limit, a maximum 

■ilu A 

Mach number contraint, and a positive altitude constraint. Figure A5.1.2 shows 

the form of this constraint. The maximum bank angle constraint a is a maximum 

m 

normal load constraint: 


Wg 

normal load in g's = ~ ~ tana 

g 

assumed to be 4. Thus a_ = 76°. 


The auxiliary equations that relate all the variables present in the equa- 
tions of motion (1.1) to (1.4) are 


(1.18) h = E-|- 


(1.19) a = 

a 


by using the assumptions cosy - 1 and y-0. 


(1.20) L = aL^ 

(1.21) D = Dq + D|_sec^a 

( 1 . 22 ) 

(1.23) ^ Cl qS 

a 

(1.24) Dq 6 Cp^qS 

(1.25) D. ' ^ 

a 


and the following constants: 
S = 49.239 m^ 

(1.26) 

W = 1.5569 X 10^ N 
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MAXIMUM 


Figure A5.1.2.-V Constraint as a 



h^O CONSTRAINT 





The example aircraft being used in this work is an early version of the F-4. 
Sonic speed c(h), density p(h), lift coefficient slope C. (M), zero lift drag 

coefficient Cp, (M), aerodynamic efficiency factor ri(M) and maximum thrust 

are tabulated functions. The numerical treatment of this data so as 

niaX 

to provide continuous values by construction of spline fits was described in 
Chapter 2 and so will not be discussed further here. 


A5.1.2 First-Order Necessary Conditions 

For this problem the Hamiltonian may be written as 


where 

P2 

Po 

(1.28) ^ 

^3 

^3 

^4 
P4 

Euler-Laqrange Equations 
(1.29) At 


A V 

= w 

— 1 

1 

o 

0 

■ DlSbc^o) + Xg9 : 

+ PiT{T 

- T ) 
max' 

+ ^2^ L -“s) + V3(V-V, 

a 

= 0 

if 

0<^<T^ax 

> 0 

if 


= 0 

if 

|a| <as 

> 0 

if 

|o| = 03 

= 0 

if 

V< V 

max 

> 0 

if 

V= V 

max 

= 0 

if 

1^1 

' ' m 

> 0 

if 

II 


9H 

3E 


w \8E 


. 2 V j ®^niax W seca 
W Y sT— ^ ^2-71- w ^ ^' 33 ^ 

' a 
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(1.30) Ap -gp 

= A V sin3 - A V cos3 
^ y 

(1.31) 5.^'=-|tt=0 

(1.32) \ = -|y=o 

Hence both A and A are constant. 

X j 

Optimality Cond i tions 


(1.33) If = 0 


(1.34) 


(1.35) 


A^V 

W 


" (2T - T^^^) 


M 

9a 


= 0 


= sec 


2 r 9^B ^D^VA^tana „ ^ J 

^ V W ^2 “T 2^4 tana 

*- a J 


9V 


Ar o 3Dp| 9Dj « fanrr 

0 = -ji (T - Dq - DLSec^a - V ^ - V ^ sec^a) - A^g + A^cosB 


3T 


9L 

, . ^ • o T " ‘max W seca a . 
AySina - “ ^2 ■ 2 9V ^3 


V‘ 


Transversal ity Condition 
(1.36) H(t^) = 0 

Also since the Hamiltonian is implicit in the independent variable, time, 
we have 


(1.37) H(t) = H(tf) = 0 
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Equations (1.1) through (1.8), (1.28) through (1.35) and (1.37) are the 
first-order necessary conditions for a minimum-time turn. 


i 

/ 

/ 

/ 


nH 

A5 . 1 . 3 Legendre-Clebsch Necessary Conditions for a Minimum (2— Order Conditions) 

The Legendre-Clebsch necessary condition for a minimum is 





«aT 

Hvt 

(1.38) 

H = 
uu 



Hva 



Htv 

«aV 

«VV 


if det H^^ = 0, then the time optimal trajectory is a singular arc. 
The components of have the following form:- 


(1.39) = 2y^ 

V = "Ta = 0 


Ar- 

(1.41) H.,t = = -# - u. 


'VT "TV V ""1 3V 


(1.42) = sec a 

00 


r 2 Dj_VA^ 

“W 


, W cos o . o 
+ P2 i ^ 

a 


(1.43) = sec-^a 


3D. 

gXg 2D^ Altana ^ VA^tano 3L^ 3 ,.^^ 

2 y _ 


r 


w 


a J 


2A g tana A^ 
(1.44) = -^-3 ^ 


2 + 2 sec^^o + V 


3V 


9V 


9V‘ 


3^0 2 

+ V sec a 
3V^ 


9^T 

T max , W seco 

- ::z— + >^2 — r- 


9V‘ 


a 


_L| 

^3lf 

9^L 

a 

a 

L 

L ot 

^3V ) 

3V^ . 


A5.1.4 Existence of a Maximum Mach Number Straight Cruise Arc 

Parsons (1972) gives a full justification for the existence of a singular 
arc with straight flight using intermediate thrust and maximum Mach number using 
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equations (1.29)' through (1.35) and (1.38) through (1.44). We will only pre- 
sent the characteristics of this cruise arc here. 

The cruise energy and velocity are given by the relation 

3V 

(1.45) = 0 


Since the maximum velocity constraint is of the form shown in Figure A5.1.2, 
the energy at which (1.45) is satisfied is unspecified at or above 29.949 Km* 

As it is expected that initial and final energies will be below the minimum 
cruise energy, this value is chosen as the cruise energy. Hence we have 


(1.46) 


= min[E] 


- 2. 9949 X 10^ m 


3V 


max 


9E 


= 0 


(1.47) V* = 




= 590.2 ms 


-1 


3V 


9E 


max _ Q 


(1.48) 

(1.49) 


a 

J* 


* = 




(1.50) = 1.2192 X 10^ m 


Equations (1.46) to (1.50) characterize the cruise arc and so at present 
we only consider a three-dimensional turn to a point or onto a line which is 
far enough away from the starting point so that a cruise arc is reached during 
the flight. The length of the cruise arc is determined so that the horizontal 
plane boundary conditions are satisfied. Thus, the initial turn to cruise arc 
and the final turn from the cruise arc can be determined separately. The solu- 
tion therefore consists of three parts: (1) initial turn to cruise arc, (2) 

flight along cruise arc, (3) final turn to required point from cruise arc. 


A5.1.5 Formulation of the Initial Turn 

As the coordinate system can be rotated so that the cruise heading, 6^ = 0°, 
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we only consider an initial turn onto the line y=0 within the cruise arc as in 
Figure A5.1.3. 

In order to facilitate the formulation of the initial turn we assume that 
the turn contains a constant radius turning segment with T* = T E = and 

/s lUaX C 

= as in Figure A5. 1*4. 

The cost function, J, for minimum time to x=0 is 


(1.51) 




V V 

max max 


V f 

time in cruise 



time in time in 
turn at turn at 


V = 


V 

max 


V< V 


max 


subject to the equations of motion (1.1) to (1.4), and the boundary conditions 
are now 


(1.52) 

E(to) = Eq 

E(t^) = E, 

(1.53) 

6(to) = 6o<0 

3(t^) = Bj* fi^ee parameter 

(1.54) 

x(to) = Xq 

x(ti) = free (strictly, 



^^^1^ ^ ^ without loss 

of generality x.p = 0.) 

(1.55) 

y(to) = >0 

y(ti) = R|^.^(l-cosB^), free 


(1.54) and (1.55) are obtained from the geometry of the problem. 


AS- 1.6 Determination of the Optimal Initial Turn Trajectory 

The variational Hamiltonian is the same as equation (1.27) and first-order 
necessary conditions and first integral as in equations (1.29) to (1,37). The 
following terminal adjoints may be obtained from (1.51): 

(1.56) X^(tj) = unknown 
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Figure A5.1 .3. -Break-up of Trajectory; Initial Turn, Cruise Arc, 
Final Turn. 



Figure /\5. 1 .4. -Constant Radius Turning Segment Within Initial Turn. 
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(1.57) 

“Rm-; n 

Ag(ti) = 


'^max 

(1.58) 

Ajt,) = ^ 

^ ^ W 


max 

(1.59) 

o 

II 

> 


With the results (1.58) and (1.59), equation (1.30) becomes 


(1.60) 


-V sing 
V 

max 


It is intended to integrate the state and adjoint equations backwards from 
time tj, hence determination of X^(tj) together with the values of the controls 
at tj^ will parameterize the whole problem in terms of 3^^. The value of is 
also in terms of 


(1.61) 




V = V 
E= E. 


max 


Controls at Time t j^ 

On the constant radius segment, bank angle for constant energy will be 


(1.62) 


a = tan 
c 


-1 


Pmax-Do-^ 


'’l 



J 


= 59 . 2 ^ 


since E = 0 


(1.63) R 


max 


mm 


v = v, 

E= E, 


max 


V 


max 


g tana 


= 2. 1195 X 10^ m 
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(1.64) 


T = "max(Wy = 10.45x10% 
(1.65) V = = 590.2 ms‘^ 


Determination of A ^(t^) 

Using (1.34) when |a| <a^ and \o\ < so that = 0 and 

unconstrained bank angle, 


WgX, 


(1.66) Ap = — « 

t nr\ 


2D|^V tano^ 


for continuity, 
(1.67) tana,, 


= tana. 


(1.68) V 


- = V 

t^ max 


Using (1.57), (1.63), (1.66) to (1.68) we have 


(1.69) A^(tj) = 


-W(l - cosBj) 
2D, V tan^a 

L. TDdX C 


Thus the initial turn problem has only the single parameter 
to satisfy a given Eq at 3 q or vice versa. 

Choice of Optimal Controls Along Trajectory 
The optimal bank angle is obtained from 

(1.70) a* = min 

where a^ is given by equation (1.17) and by equation (1.16). 
optimal bank angle, a^, is obtained from: 

(1.71) = tan 


WgA 


3 




a = a^ = optimal 


3 2 to be selected 


The unconstrained 
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Equations (1.28) and (1.33) require that: 


(1.72) 

(1.73) 

(1.74) 


> 0 if T=0 since v, = i,.|. 


A^V 


a 0 


max 


Xj. = 0 if 0<T<T 
E max 


^ Q ^ = ^max ^1 " WT 


-A^V 


> 0 


max 


Hence an observation of the sign of (determines the optimal thrust. A 
more complete (derivation of a* and T* appears in Parsons (1972). 

The optimal velocity could be obtained from the remaining optimality con- 
dition equation (1.35) with equation (1.28). This requires a one-dimensional 

search in V at constant state from the V ^ constraint to a value around the 

max 

stall velocity. However, to eliminate the need to calculate the partial deriva- 
tives of equation (1.35) in the search, the optimum velocity was obtained by 
minimizing the variable portion of the Hamiltonian: 


(1.75) V* = arg min 


A^W 


(T* “ Dq - D|^sec 


, ^ q tana* 
V 


V cosa* 




max 


constant 

state 


where V|^j^ 


is the lower limit of the search. 


A| - Equations 

Using — and A = 0 and equations (1.33) to (1.35) to get the y(t), 

Aw y 

max 

equation (1.29) for A^ becomes: 


AcW/BD. 3D, . 3T^^^ 

(1.76) " T (■^ W~ ^ " ^3E 


when la <o^ and V < 

' s max 
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(1-77) 


when |o| < 

(1.78) 
when |a| = 

(1.79) 


when I a 


3T 


_ ^ ^ 2 

■ W ^3E 3E o - Cgg 


max 


3V 


max 


3E 


cos 6 
V 

max 
3D, 


-^[e(T + v^)-(Do+v^) 


A„q tanal 




j 


°s 


Xj,V/3Dp, 3D, ^ 3T ^ 

? _ t / 0 . L c ni9X 

Af I I I r f\ r- sec Q “ ^ 


t W \ 3E 3E 


2nd ¥ < 

s m max 


ArV / 3D« 3D, 

A = -L_( __0 + __L 
'^E W ^3E 3E 


3E 


/ \WL cos' 

• 'a 


9\ 


VL sinacosa 
cos (7 a 


) 


9 3T^^, \ 3V^^^ 

sec^o - g . I + 


3E I 3E 

3D, 


CO S3 
V 

V max 
3D. 


f / \ / \ / \ o 1 ^o9 tano 

- -# [5(t+V ^)- (Dq + V ^)- (Dl + V ^)sec2c] + 


3L 


a 


3V 


\WL cos^a 
a 


9\ 


VL sinacosa 
a 



3L 

> + “ 


gAg ’ 

h 

^ 3E 

2 

WL^cos a 
L a 

VL sinacosa 
a 


a and V = V 
s max 
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(1.80) V = 


A^V/8Do _ 2 


dJ 


E W \3E 9E 


sec a “ c 


max 


8E 


9L 


a 


8E 


8L 

c 

8V 


a 


COS 3 
V 

max 


4R 


J w ^^maxN /p J. w 


) 


/ ^ *^L\ 2 


A_g tana 


r 


when I a| = cr = a and V < V and where C = 1 if T = T and ^ = 0 if T = 0. 

This completes the formulation of the initial turn problem with the single 
parameter 3j- A choice of 3-^ specifies all the initial conditions for the back- 
ward integration of the E, 3, x, y> A^ and A^ equations to produce an optimal 
trajectory to satisfy a given Eq at 3 q or vice versa. 


A5.1.7 Alternate Single Parameter Formulation 

The formulation of Section A5.1.6 shows that a = a =59.2° at t=ti. The 

c 1 

optimum bank angle is almost zero until t is very close to t, , thus producing an 

J- 

extremely sharp spike in a. This is very pronounced for |3^| < 1° x 10~ . As in 
Figure A5.1.5, since digital computations are being used, a is held constant 
over each integration step, thus broadening this spike in bank angle. To 
attempt to follow the sharp change in a more accurately, the step length must 
be reduced considerably, resulting in increased computation time. Thus, to 
avoid these numerical difficulties, an alternate single parameter formulation 
may be obtained using A^{1) as the parameter and integrating backwards from t=l 
instead of t=t^. A review of the results for 0<t<l for small |32^| revealed 
that Ar, E and x are all independent of A^ and 3. Hence the values of Ar(l), 

t ^ p t 

E(l), and x(l) for |3J =10" can be used as initial values to begin integration 

^ "6 

at t = 1 for trajectories given by |3]^| < 10" . These values are: 

(1.81) X^d) = -1.52 X 10'^ ra'^s 

(1.82) E(l) = 2.9907 x lo'* m 
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(1.83) x(l) = -0.0278 m 

The controls at t=l can be taken as 

(1.84) V*(l) = 

(1.85) o*(l) = 0 

(1.86) T*(l) = 

Further, a study of the solution near t=t^ from Parsons (1972) gives an 
approximate value for 3(1) as: 

(1.87) 3(1) = 0.53 X^(l) 

Thus we now have the initial turn problem parameterized by A^(l) to be 

selected to satisfy Eq at Bq or vice versa. 


A5.1.8 Numerical Results for Initial Arc 

Several, optimal trajectories were obtained by variation of the parameter 3^ 
or Xo(l). Parsons (1972) shows that the complete set of initial conditions is 

p 

covered by the range of parameter values 0 < 3, < •'2. 15°. For values- of (bJ above 

-3 ^ ^ 

about 10 , the entire initial turn is on the maximum velocity constraint. For 

_6 

|3i| < lO” the Ag(l) formulation was used and we present here the characteristics 

^ -17 

of the optimal trajectory obtained for a choice of Ao(l) = -10 which displays 

p 

the zoom dives that reflect the energy-state approximation used in the solution. 

Figure A5.1.6 shows the flowchart for obtaining an optimal trajectory back- 
wards in time from the cruise arc. Box A and the dotted feedback loop indicate 
logic to choose the correct value of the parameter 3^^ or Xg(l) to arrive at 
specified initial conditions Eq and 3 q. This logic would interpolate between 
the family of Bj or Xg(l) trajectories (flooding method). This may involve a 

few iterations on those parts of the trajectory where 3 changes. 

-17 

For the example trajectory chosen (Ao(l) = -10 ), the stopping criterion 

p 

used was to test if the |a| locus had been reached. This is the point 

where the angle-of -attack, a, and the bank angle are at their maximum values 
simultaneously. Figures A5. 1. 11 - A5, 1.23 give the time- histories for all the 
states and controls and for some of the derived variables.* The changes in 
velocity and altitude appear in Figures A5.1.13 and A5.1.14 where the zoom dives 
clearly indicate the instantaneous tradeoff between height and velocity. Figure 
A5.1.15 reveals that the whole initial turn consisting of 104° is completed 
within the first minute of the trajectory. The change in a (Figure A5.1.16) is 
also completed in the first minute corresponding to the change in heading angle. 

The spike in a introduced by the constant velocity turn upon reaching the cruise 
arc may be ignored in practice with no effect upon the optimal trajectory. Note 
that the entire initial turn is accomplished using T* = T^^^: the changes appear- 

ing in Figure A5.1.17 are due entirely to the changes in height and Mach number 
upon which T^^^^ depends. As expected from Figure A5.1.4 where the initial turn 
was formulated, the necessary adjustment in the y coordinate is completed in the 
first minute (Figure A5.1.19) during the change in heading angle. Figure A5.1.21 
shows the time-history of the derived variable a. Notice that a=12° at the 
beginning of the turn where \o\ " which is the max turn locus. Figure A5. 1.2-1 

*The Figures A5. 1.11 - A5. 1.42 are collected at the end of Appendix 5.1. 
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Figure A5.1 .6. -Flowchart for Optimal Trajectory. 








is an altitude-range profile of the trajectory where the very fast. changes in 
altitude are clearly visible. The variations in Mach number with range appear 
in Figure A5.1.25. The various combinations of initial conditions (3 q»Eq) that 
can be reached with this ^^(1) value can be seen in Figure A5.1.26 while Figure 
A5.1.27 shows changes in a with energy-height. The horizontal plane projection 
of the initial turn trajectory appears in Figure A5.1.28. 

A very interesting view of the trajectory is obtained from the altitude- 
Mach number profile of Figure A5.1.29. The energy contours have been added in 
as well as the |cr| maximum velocity constraint and the a=0 with stall 

constraint. Starting on the maximum turn locus, the trajectory consists of a 
dive at approximately constant energy till the velocity has increased to about 
0.9 M. The velocity remains approximately constant at 0.9 M while the trajectory 
gains altitude. Using Figure A5.1.26 which gives the change in 3 with energy 
together with the altitude-Mach number profile, the progress of the trajectory 
can be interpreted with regard to the change in heading angle as well. By the 
time the heading angle relative to cruise is below about 5°, the trajectory 
essentially coincides with the A3=0 min-time path of Parsons, Bryson and 
Hoffman (1975). The trajectory follows this path to the constraint and 
then at V to the cruise point. 


A5 . 1 . 9 Formulation of the Final Turn 

There are two general types of final turns: (1) those beginning with a 

period of constant energy-maximum thrust turning and (2) those beginning with 
minimum thrust and straight bank angle chatter. Each of these is now described 
in detail . 

(1) Turns beginning with maximum thrust . The approach here is similar to 
that for the initial turns. We assume an initial turning segment with T=T , 

^ J I IQ A 

E=E and V=V ^ as in Figure A5.1.7. 

C iTlaX 

The minimum time cost function may be expressed as: 


min 

t,a,V 





max 


max 


V 


f^f 

Jt„ 


dt 


time in cruise 


max 


time in 
turp at 
V= V 

max 


time in 
turn at 
V< V 

max 
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subject to equations of motion (1.1) - (1.4), (1.12) - (1.15) and the boundary 
conditions 


E(t2) = 

E(tf) = Ef 

6(t2) = 62* parameter 

6 (tf) = 


X(tf) = Xf 

y(t2) = R„^n(l-cos62). f''ee 

y(tf) = yf 


The Hamiltonian for the problem is as in equation (1.27). The first order 
necessary conditions and first integral are given by (1.29) - (1.37) . We now 
have the following initial ad joints derived as for the initial turn: 


\i^2) 

V^2> 


-w(l-cos$2) 

- COS62) 


max 
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A y 

max 

The determination of the optimal final turn trajectory closely parallels that 
of the initial turn as given by equations (1.60) - (1.80), with appropriate sub- 
script changes. The state and adjoint equations are integrated forwards from 
time t2 and the problem is parameterized in terms of 62* At each step the opti 
mum controls are evaluated using equations (1.70) - (1.76). 

(2) Turns beginning with minimum thrust and bank angle chatter . These 
turns require an initial period of straight flight with minimum thrust and bank 
angle chattering between positive and negative minimum values. This is a maxi- 
mum deceleration arc and is similar to a singular arc as it only occurs with 
Ag=0. The Hamiltonian with A^="0 is given by: 

X.V 2 

H ** 1 + (T - Dq - 02Sec a) + cos6 + X^V sin$ 

which is independent of the sign of the bank angle. Application of the Minimum 
Principle yields: 

T* = 0 (X^ > 0) 


= *°chat<V) 


V* = arg minfH = 1 - cos6^ + X^V 


'^chat^'') ' min[a5(V).0j 
= cruise heading 

Note that as a result of our modeling, bank angle can be changed instantaneously 
and so an equal time chatter between positive and negative bank angle can be 
made so that 6 = 0 on the average and hence X^ remains zero. For a better under- 
standing of the geometry of the problem consider the Minimum Principle applied 

to the problem with 6=0 so that X = 0: 

^ y , 
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im*n = (Xft.Ap,A^)-(3,E,x)^ 

T,a,V t X 

For a given state the state-velocity space (hodograph) appears approximately as 
in Figure A5.1.8 considering just a single section in each plane. This approach 
was. an, aid in the ^determination of optimal velocity when T*-0. As In Figure . . 

A5. 1.8, the admissibly state-velocity space has a non-convex nature so that after 
leaving the. gruise arg. at. V*=. the solution can be expected to jump at some 

point along the trajectory from thy T*=0, V*=V„^„, a* =a corner to the T* = 0» 

' .max \ ■ m V-:- ■ 

V* for \a\ 1^*1 -^s“% of the state- velocity space. This know- 

ledge, reduced the onerdiwnsional . search in velocity described in equation (1.75) 
to. just a check at these two velocities. As a result computation time was greatly 
reduped. . 

The turns with bank angle chatter may be formulated as follows: 




I 





\s 

V 

\ 

s 

\ m 

n ^ y 



; o' 

N 

T=0 

v/ ^"^max 





Figure A5.1 .8. -Sections Through State-Velocity Space (Hodograph) . 


with the cost function expressed as: 


min \ 

fo= "0 ^r^idE 


+ 

p^dtl 

T,V,0 1 

V Je E 

L max c 

V= V* 


'h J 


T*=0 
a*= ±0 




chat 


time in time in chatter time in 


cruise 


turn 


which is 


min ^ J = 
T,V,0 


1 - 

E, V V 


^1 ^ r 1 \_j; 

V Je e 
max c 


V 

max, 


dE 


T*= 0 
0 * = ±0 

V= V* 




+ / dt > 
1 


chat 


subject to equations (1,1) - (1.4) and (1.12) - (1.15) and the boundary conditions 
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E(tf) = Ef 


E(ti) 

B(t^) 


Ej, free parameter 
0 


x(tj) = Xj(Ej), free 
y(tj) = free 


B(V) = Qf 
x(t^) ?= 

y(tf) = 


The Hamiltonian and first-order conditions are the same as before: (1.27) and 

(1.29) - (1.37) with the initial adjoints 


A£(ti) 


(\’Vax) 


T*= 0 

o* = ±a 

V= V* 


chat 


AoCti) = 0' for continuity across juncture with chatter arc 

p i 

\{U) - ^ 

X i V 

max 

Ay(t^) =0 


Again we have the problem parameterized in terms of a single parameter 


A5.1.10 Numerical Results 

Parsons (1972) presents numerical results using the theory described above. 
Figure A5.1.9 gives the family of trajectories in the beta-energy plane and the 
regions marked I to VIII are described by means of the control programs in 
Figure A5.1.10. For the chatter cases t^ has to be determined by integration 
of the equation: 


V= V* 

T*= 0 

^* = ±^chat 



185 



CHANGE IN HEADING ^deg.] 





Figure A5.1.10 

Summary of Control Programs for Variable-Altitude Final Turns from Cruise Arc 


00 








and then the four states and two adjoints were integrated from t^ to t^. As in 
Figure A5.1.9 the zoom climb when V* jumps is followed after a short while on 
the max-turn locus by a switch to from T=0. As explained earlier, com- 
putation time here was considerably lower than that in the initial turns because 
of the knowledge of the non-convex nature of the admissible state-velocity space 
reducing the search for V* to just two points. 

Trajectories beginning with maximum thrust were of the same form as those 
with chatter except that there is no straight flight portion. In the altitude- 
Mach space, the trajectories move down the maximum velocity constraint and then 
zoom to the max-turn locus as for the example case. The constant energy turning 
is with T ” "^o^^owed by T = 0 on the max-velocity constraint through the zoom 
and onto the locus. After a short stay on the locus, thrust switches back to 
max. These turns are suitable when the final energy is high, that is, around 
24 - 28 Km. 

A5.1.11 Computational Efficiency 

In obtaining the example trajectory of Figures A5.1.11 to A5.1.29, a ratio 
of about 2:1 resulted for CPU timeireal time. It is possible to reduce this 
ratio to at least 1:2 by making the program more efficient. It is worth men- 
tioning that substantial improvements were made to arrive at the value 2:1 as 
initially the ratio was of the order of 10:1. This was achieved by allowing a 
tolerance of 15.2 m and 0.001 M on the spline fits for all the tabulated data, 
i.e. , for changes in altitude and Mach number less than those above, new values 
of aerodynamic or atmospheric data were not obtained. The sacrifice in accuracy 
was minimal. Another time-saving method was to increase the storage in the com- 
puter (e.g.-, spline coefficients) so as to avoid calling standard spline fit sub 
routines too many times. 

In analyzing the distribution of CPU time among the major computations re- 
quired, it was found that the calculation of V* was heavy on time. Hence, the 
frequency of this computation was reduced to every 4— step which resulted in a 
considerable decrease in CPU time. The sacrifice in accuracy began to show up 

j. L 

when the frequency was reduced to less that every 4— step. In order to under- 
stand the structural form of the Hamiltonian F(v), the curve was plotted out at 
various stages along the trajectory. Figures A5. 1.30 - A5. 1.42 cover the major 
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changes in F(V) and hence V*. As expected, these occur along the zoom maneuvers 
where velocity is traded for altitude at constant energy and in very little time. 
An anticipation of when these jumps in M* occur (i.e., when zoom maneuvers are 
made) will allow a reduction in the range of V over which the search is made 
for V* at points along the trajectory between zoom dives. 

For real-time implementation (Subsection 5.5.2) computation time is con- 
siderably reduced by storing parts of the trajectory (e.g. , the A3=0 minimum 
time path). 
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Figure A5.1.16 Bank angle time history 
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Figure A5.1.21 Changes in alpha with time 




150.00 200.00 

TIME (sec) 


250 


Figure AS. 1.23 A time-history 
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Figure A5.1.41 F(V) along AB = 0 path 






APPENDIX 5.2 

HIGHER ORDER APPROXIMATION FOR h-LAYER SOLUTION 

The exact solution of the a,y controls for the h-boundary layer requires 
solving three simultaneous, nonlinear equations for a,Y and a Lagfangian multi- 
plier A. These equations are 

(2.1) 0 = (L^a + T sina) - mg cosy 
(respectively equation (5.6.2) of Section 5.6) 

(2.2) 0 = + Amg^ " ^E mg .cosa - Dq - a^nL^) + A(L^a + T sina)]cosy 

(respectively (5.6.6)) 

(2.3) 0 = -AS ~ (T sina+2riL a) + A(L +T cosa) 

^ E mg a ' 'a ' 

(respectively (5.6.7)). 

Note that the validity of equation (2.3) depends on the optimal a not occurring 
on one of its constraints (i.e. , 0 < a < must be true). 

We assume a power series expansion for cosy and A which solves (2.2), (2.3) 
simultaneously. That is, suppose that if cosy=c(a) and A=A(a) solve (2.2), 

(2.3) , then 

2 I' ‘ 

(2.4) ' c(cx) ~ Cq + c^a + + — 

(2.5) - A(d) “ + A2a + .“1. 

We will solve for c^., A^ up to i = 2. 

For i = 0, equations (2.2) and (2.3) become 

(2.6) 0 = Vs) - ^ (T-D(,)]co 
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(2.7) 0=Xp(L„+T) 

Hence, we have 

(2.8) Aq = I ' 

X 

V 

(2.9) c„ = f- 

For 1 = 1, equations (2.2) and (2.3) become 

(2.10) 0 = A^mg - [1 + — (T- Dq)]Cq 

0 = -^1 * 2nL„) +Xi(L„ + T) 


Hence, we have 


( 2 . 12 ) 

(2.13) 


A, = 


c, = 


(T+2nL ) 

E mg g 

(L +T) 

' ot 

A|V(T+2tiL^) 


For i = 2, equations (2.2) and (2.3) become 

(2.14) 0 = Xging * [1 + X| (T- 0 q)]c2 - [-X| (7 + hL„). + 

(2.15) 0 = X2(L^+T) 

Hence, we have 


(2.16) A2 = 0 

Equation (2.14) becomes after substituting A^ from (2.12) 


^l(l.„+T)]Co 

I 
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(2.17) 0 = [1 + 




S?'T-Do)]c2 


E mg 'Z 




Hence, we have 


(2.18) 


C2 


yvV*X(T ,,L ) 

\V^y E mg '2 'a' 


The preceding perturbation series for y(ct) is used to solve for a from (2.1) 
as follows. Substitution of (2.4) in (2.1) gives 


(2.19) 0 = L^oi + T sina - mg(cQ + C|«+ C 20 t^) 

To be consistent, we should expand sina to second order in (2.19) and solve 

(2.20) 0 = L^a + Ta - mg(cQ + c^a + C 2 a^) . .. 

We could solve (2.2) directly by quadratic formula. Instead we will expand a in 
a perturbation series with respect to the parameter 

(2.21) e = 

a 

Thus, (2.20) becomes 

(2.22) 0 = a - e(Cq + c^a + C 2 cx^) 
where we assume a(e) is given as ' ' 


. ' ' r - . ri: 

(2.23) a = Uq + ea^ + e^a2 + ... 

For i = 0, (2.22) becomes 

(2.24) 0 = ciQ. .. , . _ . . ' , . .. ^ 
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which gives ag directly. 

For i = 1, (2.22) becomes 

(2.25) 0 = - Cq 

which gives as 

(2.26) aj = Cq 

For i=’2, (2.22) becomes 

(2.27) 0 = 02 - (c^oi^) 

which gives ' 

(2.28) 02 = CjCq 

Thus, to second order in o we obtain 

(2.29) o - CqC + c^CqC^ 

In terms of the aerodynamic coefficients (2.29) gives o as 


w- mg 


(2.30) 


o - 


j- (mg)^|V(T + 2nLj 
c 


^ (T - T] [1 + (T - Do)]2[L„ + T]3 


To obtain y we may solve (2.1) directly for y, using the o from (2.30) or we 
may use the series expansion for cosy in (2.4) and solve for y from that, ^ Ig- 
noring terms with 0(e) dependence would give 


(2.31) 


y - cos 


-1 


L + T 

o 

mg 


o 


where o in (2.31) is given by (2.30). 

Note that at this level there is no reason to work with y rather than cosy, 
and one may as well avoid cos”^. 



The two things we need from the h-boundary layer are y* (or cosy*) and AJj. 
We have obtained the former in (2.31). We now obtain Aj!j. The adjoint Aj![ is 
such that 

(2.32) , 0=1- COSY* + A| (T .cosa* - Dg - nL„(a*)^) + Af;V siny* 

(respectively equation (5.6.9)). 

By rewriting equation (2.2) we obtain 

(2.33) {1 - ^.cosy + (T cosa “ Dq “ a^nL^)}cosy 

2 V 

= (1-cos y)(-^+ Amg) + A(mg cosy - L^a - T sina)cosy 


Using equation (2.1) we can simplify (2.33) to obtain 

(2.34) (1 - ^ COSY + ina J}cosy = s1n^y(^ + Amg) 

c ^ c 

comparing (2.34) to (2.32) we see that A and AjJ are related as follows: 

(2.35) A|J = - (X + A*mg) 

c 

where A* is obtained from (2.5) with a* substituted from (2.30). Note that A* 
is given as 

(2.36) A* = A^a* 


and from (2.12) we have 


(2.37) A* = 


+ 2t)L ) 

E mg ^ 

L + T 
a 


a’ 


and hence Ajj is given by 
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(2/38) 


X*= - 


tanv*/V , . 

: , (L + T) ■ 

a 


Thus, we may use either (2.38) or (2,32) directly to obtain XJj from y* and a*. 
The direct computation of XjJJ from (2.32) may be more accurate but (2.38) seems 
to give a better idea of how behaves in terms of V, y*, a*. Note that the 
only error in (2.38) comes from the approximation of X*, and this approximation 

q 

has an error of 0([a*] ). 

Constraints . The next possibility to check is whether or not a or y takes 
a value on a constraint. To denote possible constraints, define y^ as the posi- 
tive solution Y of 


(2.39) 0 = + T sina^ - mg cosy 


if such a solution exists, i.e. , if we have 


L + T sina. 

(2.40) ^ 1 

mg 


Likewise, define a as the solution a of 
m 


(2.41) 0 = L^a + T sina - mg 


Analysis of the original optimization problem, to minimize (for y > 0) the 
function 


1 - cosy 


(2.42) 


X*-^ 

mg 


(T 


V siny 

(respectively, equation (5.6.4)), 


2 

cosa - Dq - nL^a ) 


subject to the constraints 


(2.43) 0<Y<7T 

(2.44) 0 < a ^ a^ 

and (2.1), shows that as y-^-O, the function in (2.42) tends to +~ as Hence, 
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we must have 


(2.45) a < a 

' m 

where is defined by (2.40) above. Note that ci^<a^ or niay occur. 

Thus, the only constrained values of (a,y) are due to the constraints on a 
in (2.44). 

If then we check the possibility of 

V V / 2^ 

1 - cosy. + AS — (T cosa_ - D« - nL a_ ) 

V 's Emg s 0 'as' 

(2.46) -A? = £ 

V siny^ 


being optimal. Likewise, we must check the possibility of 


n 1 + ^ (T- Dn) 

(2.47) -X? = ^-523 — 

V 


being optimal (the a=0, y=7t/ 2 solution). 

c 

The optimality test for these constraints simply involves comparing Au and 

0 ^ 

Ah to A^ for the unconstrained problem. The comparison goes as follows: 

Case 1 . a*<0 or a^<a*, either A^ or A^ is optimal. The optimal case corres 

ponds to the larger or the two (remember -A? or -A? is minimum). 

Os n n 

Case 2 . 0<a*<ag, either A^, A^ or Ajij is optimal. The optimal case corres- 

ponds to the largest of the three. 

h-Boundary Layer Algorithm 

1. Calculate from (2.41). 

2. Calculate a* from (2.30) ( second order approximation ) or from the first 
order approximation 


V 

w- mg 


(2.48) 


- 


[1 + 


yk JL. 

mg 


(T-Dq)][L^ + T] 


or, if you prefer, by solving the nonlinear equation 
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(2.49) 


0 = L a + T sinot ~ 
a 


V 

y-mg 

c 


1 + A* ^ (T - D.) 
t mg u 


by a Newton-Raphson method. 

3. If then go to 6, 

ms’ ^ 

4. If a* < 0 or a < a* go to 7. 

0 ^ s 

5. Calculate A|^, X^: the largest corresponds to the optimal solution, 

X^ a = 0, y - u /2 
Xjj ^ a*, y* 

Go to 8. 

6. Calculate a*, if a* < 0 then X? is optimal and a=0, y = tt/2. If a*>0 then 

0 ^ 

calculate A^, AJj and choose the larger •'-this gives the optimal. 

Go to 8. 

0 s 

7. Calculate A^, A^. The larger corresponds to the optimal solution. Go to 8. 

8. If h<h(E)> keep y>0 as it is (y = Y 3 > 'fr/2 or y*). If h>h(E), take -y 
instead. 

STOP 
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APPEND(a s.i 

HIGHER ORDER APPROXIMAnON OF ^-LAYER SOLUTION 

This appendix presents a hiahei appi oximation of the y-boundary layer 

calculation given in Section ‘■..7. Tl)e in-incipal improvement is to use a Taylor 
series expansion to determine pprt'jrbntion of the angle -of -attack from the 
value which maintains a steady flinht patP male. By calculating 6a, the per- 
turbation of the anqle-of-attacL from its steady state value instead of 
calculating the angle-of-attack i directlv, we vn‘11 obtain better numerical 
behavior of the solution. 

To define the problem more precisely, we must make some preliminary defi- 
nitions. Let L(a) be the function of the anal e -of -attack o. defined by 

(3.1) L(«) = 1 - cosy + * n cfji.. - DQ-n^aR + X*V siny 

c 

and let f{a) be the function of defiii»^d hy 

(3.2) f(a) ^ L^^^o: + T '^ino - mg cns> 

Note that the functions L end f alco depend on the flight path angle y- The 
optimization problem is to mi ui-i't -e r'\o rak io L/f if y '' ,*{E,h) and to maximize 
this ratio if > ,*(E,h). wLei - 'U -leiiotes the pseudocontrol value of y 

calculated in the h-boundarv laye*- lu toe minimization v/e are to restrict 

values of a such that f(- ) 0; lUev.'ise, in the maximization we are to restrict 

a SO that f (a) 0. In addition ar-^ inequality constraints on a, namely 

(3.3) 0 < a < 

Let us define the steady state oT as rfie value of a, denoted a^, which 

solves the equation 

(3.4) 0 = f(..Q) 

That is, is the anal p-‘>f att [ ..• •i ' m-, in • in steady state. Note that 
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depends on y (as well as E and h through the aerodynamic coefficients and 
T). As pointed out in Appendix 4.2 on the existence oF solutions to boundary 
layer calculations, a solution to the optimization problem will exist in this 
case with a>aQ fory^Y* and a<aQ for The former case corresponds to 

f (a) > 0 and the latter case corresponds to f (a) < 0. Thus, in either the minimi- 
zation or the maximization it will be necessary for a to be found in the interval 
ttQ<a<a^ (for the minimization) or in the interval 0<a<aQ (for the maximiza- 
tion). In each case, the optimal a is either an unconstrained value or it is 
(in the case of the minimization) or 0 (in the case of the maximization). 

As we pointed out in Section 5.7, when E < and y>y* the optimal a is 0; 

similarly, vihen E>E^ and y<y* the optimal a is Thus, vie will consider only 
the other two cases when E > and y>Y* or when E<E^ and y<y*. To determine 
the optimal a in these cases we set the derivative of L(a)/f(a) equal to 0 in 
order to determine the unconstrained values of a and then compare these values 
to the corresponding constrained values of a. Setting the derivative of L(a)/f(a) 
equal to 0 gives the equation 

(3.5) f(a)L‘(a) - f'(a)L(a) = 0 

r 

To obtain an approximate expression for the solution (unconstrained a) of 

(3.5) we expand L(a) and f(ct) in a Taylor series around a^. Thus, we obtain 

(3.6) L(a) = L(a^) + L'(aQ)6a + J5 L"(o:q) (6a^) 
and 

(3.7) f(a) = f’(aQ)6a + («q) (<Sot^) 

Note thac we have neglected terms higher than second order. Substituting (3.6) 
and (3.7) into (3.5) and collecting terms gives the follovnng quadratic equation 
for a. 

(3.8) f (ao)L(a(j) + f'(ag)L(ao){6a) + !s(f"(aQ)L' ( oq) - f ' (aQ)L” (ocg)) (6a^) = 0 
The solution of this equation is given by the quadratic formula 
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(3.9) 6a = 

-f"(ag)L(aQ) ±-V^l(f"(“o)L(ao))^-21''(aQ)L(ag)(f”(aQ)L’(aQ)-f’ (ag)L"(ctg))J 

{f‘(ag)L’(c<o) - r(aQ)r(ag)} 

To determine which root to choose we must calculate the derivatives of f and L 
from the original expressions (3.1) and (3.2) for f and L. Thus, we obtain 

(3.10) L(qq) = 1 ~ COSY + j~ (T coSclq - ^Q-nL^aQ^) + AjiJV siny 

(3.11) L'(ag) = -A| ^ (T sinag + anL^Cg) 

(3.12) r(ag) = (T cosag + 2r,g 

(3.13) f(og) = + T cosag 

(3.14) f”(aQ) = ~T sinaQ 

Substitution of (3. 10) -( 3. 14) into (3.9) gives us 

(3.15) 6a = - - 



v/here K is given by 

A*V 2 

(3.16) K ^ - — (T"^ + 2nL^T(aQSinaQ + cosaQ) = L^T cosa^ + 2n) 

and L is LCag). Note that K has the opposite sign of A^. Also note that L, T, 
L^, cosa^, sina^ are all positive. Knowledge of the signs of these quantities 
will enable us to determine which root to take in (3.15), 

Suppose that E < and y<y*. In this case we must have 6a >0, Since 
E<Ec> the adjoint is negative and thus the coefficient K is positive. Thus, 
in this case the unconstrained value for 6a is given by 


Isin^a, 


- 2(L + T COSar,)^ 
a 0 TL 
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On the other hand, if E > and y>y*, then we must have 6a <0. Moreover, the 
coefficient K is negative in this case. Unfortunately, this means that both 
roots of (3.15) are negative and we cannot immediately distinguish which is the 
correct root. Moreover, for K negative it is possible to obtain complex roots 
in (3.15) and we must be able to handle this case adequately to obtain a real 
value for 6a. What follows is an algorithm for calculating 6a in the various 
cases of E and y- ' 

Y“Boundary Layer Calculation: Algorithm 

1. If E < and let a=0 and stop. 

2. If E > E^ and y < Y*» let a = and stop. 

3. If E < and Y<y*, calculate 6a from (3.17) and go to step 5. 

4. If E > E^ and Y>y*, calculate both roots from (3.15) and go to step 6'. 

5. If aQ^6a>a^, let a = a^ and stop; otherwise, go to step 7. 

6. Replace complex roots by real part. If UQ+6a < 0 (for both roots) , let 
a=0 and stop. Otherwise discard any root such that aQ+6a<0 and go to 
step 8. 

7. Let = ciQ <5a. Choose a=a^,a^ whichever minimizes L(a)/f(a). Stop. 

8. Let a^ = aQ+6a (both roots for 6a). Choose ot = a^ or a=0 whichever maxi- 
mizes L(a)/f(a). Stop. 
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CHAPTER 6 

REAL-TIME SPT CONTROL LAW 


This chapter describes the real-time SPT control law as it has been im- 
plemented and tested on the CDC 6400 computer. Section 6.1 contains a descrip- 
tion of the control logic for different parts of the flight, namely before the 
cruise arc (Case I), on the cruise arc (Case II), and after the cruise arc 
(Case III). An estimate of the total computations required is also given. The 
exact computer time for each cycle of the control computation will depend on 
the characteristics of the flight computer. This in turn will determine the 
control update rate, which can vary for different control loops. In general, 
the computation rate for the faster layers must be kept higher than that for the 
slower layers in accordance with the SPT approximation. Fortunately, the com- 
putations involved at faster layers are less than those at slower layers so that 
the SPT algorithms of Chapter 5 are easily implemented to satisfy this require- 
ment. 

Other sections of Chapter 6 contain results of using the real-time control 
algorithm for an F-4 aircraft. The results are very encouraging both in terms 
of accuracy and computation time. 


6. 1 Description of Real-Time Control Logic 

The purpose of real-time control logic is to compensate for changes in 
target relative position, velocity and heading. A change in target velocity 
produces a change in the intercept point, namely, t^, x(t^) and y(t^). If the 
pursuer's position is on or before the cruise arc, this requires a change in 
cruise heading 3^ and time t^ for coming off the cruise arc. It is convenient 
to separate the control logic into three parts depending on the position of the 
pursuing aircraft relative to the target. 

6.1.1 Case I: Before the Cruise Arc 

If the pursuer is on the energy climb path to the cruise arc and is 
changed, a change is required in 6 before the cruise arc. For small 3- changes, 
a linear control law relating o to (E,3) can be used. This may produce small 
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changes in (E,h,y) which again can be compensated by linear feedback law. The 
aircraft state would return to the energy climb path after the required correc- 
tion in 3 has been made. 

If the change in 3 is large, Parsons' (1972) energy state solution suggests 
the use of nonlinear control law. This requires zoom climb or dive to the max- 
turn-rate (MTR) locus, heading change through the appropriate amount and zoom 
dive or climb to the min-time energy (MTE) climb path. Boundary layer control 
calculations for h and y required to complete the zoom maneuvers. The 
latter are basically vertical plane maneuvers at constant thrust, requiring cal- 
culation of a to complete the zooms. However, simulation results for the F-4 
have shown that it is better to turn at current energy without zooms since the 
zoom dynamics is quite slow (see Tables 6.2.1, 6.2.2 and 6.2.3). 

6.1.2 Case II: On the Cruise Arc 

Control law similar to Case I is used here except that it is much simpler 
in form since the linearized system is time invariant and the optimal cruise 
energy is constant. Thus, constant feedback law is appropriate for small changes 
in 3j,» h and y. 

6.1.3 Case III: After the Cruise Arc 

The terminal part of this case involves short range optimization problem 
which is beyond the scope of the present study since SPT is not valid for this 
case. However, for small changes in terminal conditions, a linearized solution 
is appropriate and this would take the same form as the linear solution in Case I 

6.1.4 Feedback Structure of the On-Line SPT Algorithm 

Figure 6.1.1 shows a schematic of the on-line SPT algorithm. From a feed- 
back control viewpoint, it is a hierarchical structure with six loops. The 
fastest loop, i.e., y-loop, is at the bottom of the hierarchy and it requires 
inputs from all the higher loops. Since the higher loops operate at slower 
speeds, these inputs to the y-loop are updated at different rates which are much 
less than the speed of the y-loop. The command y* is provided by the h-loop 
and the nonlinear SPT control law computes a to follow y*. If the error |y-y*| 
is small, a linear feedback control law involving both h and y is used to take 
into account the interaction between the two loops. 
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Target position 
and velocity 



Figure 6. 1. 1 

Hierarchical Feedback Structure of the On-line SPT Algorithm 


x-y loop 


E-loop 


6-loop 


h-loop 


T-loop 













The sampling times and the control update rates of the different loops in- 
crease as one goes down the hierarchy. This causes no problem for the SPT al- 
gorithm since the optimization problem at lower levels is simpler than the opti- 
mization problem at higher levels. In Figure 6.1.1, Case # corresponds to 
Cases I, II and III. The notation a|J denotes the optimal value of a in the h- 
boundary layer, which is the asymptotic value for a in the y boundary layer. 

Thus the calculations of a at higher levels serve as checks on the calculation 
of a at lower levels and may be used to shorten the duration of the boundary 
layers. 

6.1.5 Real-Time C apability Assessmen t 

In this section, we provide execution time and storage estimates for the 
different control loops of Figure 6.1.1 on a Texas Instruments 9900 microcomputer. 
TI9900 has a memory of 32K, 16 bit words and basic instruction execution times 
in microseconds are given in Table 6.1.1. The CPU times per iteration, nominal 
number iterations and storage requirements for control loops of Figure 6.1.1 
are given in Table 6.1.2. The Y”loop, for example, takes 1.5 msec of CPU time 
per iteration for a control a update and, therefore, can be iterated up to six 
times for 9 msec if the measurement sampling interval for y is, say 10 msec. 

Slower sampling rates for y state would allow time for more iterations of the 
control loop. The control update rate for the complete on-line calculation 
would be determined by the radar measurement rate. Assuming 10 samples per sec- 
ond, a total of 100 msec are available for control computation. If the CPU time 
is allocated equally between all the five loops, 20 msec of computation will be 
available for each control loop. This implies that two iterations of the y-loop 
calculations can be performed 6 times during 20 msec so that a sampling rate as 
high as 250/ sec can be used for y. It is clear from these figures and the num- 
bers in Table 6.1.2 that the proposed algorithm can be easily implemented on- 
line using a TI9900 microcomputer. 

The estimate of 27 msec for a complete control calculation of all the loops 
including linear feedback control implies that radar measurement rates of as 
high as 30 samples/second can be used. This would allow interception of highly 
maneuverable targets. On the other hand, if lower radar measurement rates are 
used, CPU time is available for further refinement of the SPT controls. In 
particular, the implementation of one or two iterations of the "continuation 
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Table 6.1.1 

Execution Times for TI9900 


Instruction 


Execution time 
(microseconds) 

(Clock rate is 3 MHz) 


Branch 

Register to register 


2.67 


Add (words/ bytes) 

Register to register 
Indirect to indexed 


4.67 

8.67 


Multiply 

Register to register 
Divide 

Register to register 

Shift (left/right) 

1 bit 
8 bits 

Move data (words/ bytes) 
Register to register 
Register to directory/index 


17.33 

41.33 

4.67 

9.33 

4.67 

7.33 


Load communications register unit 
(register to CRU) 

8 bits 

16 bits 

Store CRU (CRU to register) 

8 bits 
16 bits 


12 

17.33 


14.67 

20 







Table 6.1.2 

TI9900 CPU and Storage Requirements for 
the Real-Time SPT Algorithm 



Number of 
average 
iterations 
required 

CPU/ 

iteration 
w (msec) 

Average 
cycle time 
(msec) 

Storage 

(words) 

- calculation ' 

4 

2 

8 

400 

3^ - calculation 

4 

1 

4 

400 

E - BL 

4 

1 

4 

210 

Linear feedback around x- 
solution (u,cr,a controls) 

•y 

' 1 

0.5 

0.5 

20 

Linear feedback around 
E-BL (a, a controls) 

1 

0.5 

0.5 

100 

3 - BL 

2 

2 

4 

500 

Linear feedback around 
B-BL (a, a controls) 

1 

0.5 

0.5 

100 

h - BL 

4 

0.5 

2 

10 

Y - BL 

2 

1.5 

3 

10 

Linear feedback for 
h and y {a controls) 

1 

0.5 

0.5 

400 





2150 

Estimated storage for 
software 




4000 

TOTAL 



27 

6150 


BL boundary layer 

Typical sampling interval for radar returns = 100 msec 
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algorithms" discussed in Chapter 3 seems feasible within a total CPU time of 
100 msec. Table 6.1.2 also shows that the storage requirements , are well within 
the 32K limit of the TI9900. 

6.2 Simulation Results 

In this section and the next, we describe simulation results obtained using 
the real-time SPT control algorithm. First a more detailed description of the 
actual algorithm is given since certain minor but important modifications to the 
SPT control law were required. 

6.2.1 Climb to Cruise Arc 

The optimal flight path to the cruise arc is obtained by retrieving from 
storage V*{E) which provides the necessary information to construct the E boun- 
dary layer, h and y boundary layers are then added onto the E-layer solution. 

6.2. 1.1 E-Boundary Layer Calculations 

As mentioned in Subsection 5.3.2, the necessary quantities for the E boun- 
dary layer control problem are V*(E), h*(E), a*(E), u*(E) and A*(E). Having 
determined V*(E), however, the remaining quantities are directly determined from 
E and V*(E). 

V*(E) was determined off-line for a fine grid (1 point/61 m) of values of 
E, and a piecewise linear approximation was then applied. For E<E^j V*(E) was 
determined by maximizing (5.3.10) subject to (5.3.11). The maximization was 
performed by exhaustive search on a grid (1 pt/(12.2 m/sec)). For E>E^, V*(E) 
was determined by minimizing (5.3.10) subject to ( 5 . 3 . 11)3 again by exhaustive 
search. For E>E^, it was found that V = V^= 590.2 ms”^ was optimal. 

Figure 6.2.1 shows the piecewise linear approximation for V*{E). The approxi- 
mation was chosen such that its divergence from V*(E) was always less than 

6.1 m/sec. We note that there is a jump from 585.2 to 590.2 ms”^ at E = E^. 

Figure 6.2.2 contains a flowchart which represents the E boundary layer 
control program. This flowchart may be regarded as a detailed version of Fig- 
ure 5.3.1. 

dE 

We note that in the last block of the flowchart the quantity ^ is computed. 
This is the implementation of (5.3.1) with a = a*, V = V* and iJTjyjax " 
calculation will not be required on board the aircraft since E would be measured 
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V*(E) (ms 






Determine V*( 
linear approxima 



Compute aerodynamic 
coefficients given E,V*(E) 




6 . 2.2 

lary Layer Computations 
















or estimated from sensors. 

A simulation of the aircraft trajectory has been run using the controls 
calculated for the E-boundary layer. The simulation consisted of a predictor- 
corrector (Adams-Bashforth , second order) solution to equation (5.3.1), with 
a= a* and determined by the control algorithm at each step. The stepsize 

was fixed at .5 seconds, and the simulation was run past the point where the 
aircraft reached the cruise energy-height. The initial energy-height was 5482.7 m. 

Figures 6.2.3- 6.2.9 show the computed time histories of energy-height, 
altitude, velocity, Mach number, thrust, range and a respectively for the E- 
boundary layer trajectory of the aircraft. 

The aircraft, with these boundary conditions, reaches cruise energy after 
359 seconds. This is slightly slower than the 351 seconds taken by aircraft 
with the same boundary conditions (zoom changes in height and velocity are assumed 
instantaneous) in the energy state approximation of Section 6.3. This dis- 
crepancy is due to the cruise approximation of V*(E) by five piecewise linear 
functions. 

Figures 6.2.3, 6.2.4, 6.2.7 and 6.2.9 show some chattering behavior on the 
cruise arc. This corresponds to oscillation around the cruise arc which can be 
eliminated by adding some anticipatory action as the cruise energy is approached 
or by switching to a linear multivariable control lav/ in E, h and y. 

Aside from the slight differences in time-unti 1 -cruise , Figures 6.3.1 and 
6.2.3 show energy-height time histories which are quite similar. The time his- 
tories of altitude, velocity, Mach number and range are also quite similar for 
the two algorithms, except at the beginning, where the real-time E-B solution 
requires a zoom-dive to the MTR locus. The thrust histories are slightly dif- 
ferent, being smoother for the present (E-boundary layer) case. Angl e-of-attack 
(a) appear to differ between the two solutions, but the effect of a is small. 

We have also included in Figures 6.2.10-6.2.12 profiles of altitude versus 
Mach number, altitude versus range, and Mach number versus range. These are 
also similar to the corresponding plots for the energy state trajectories (Fig- 
ures 6.3.12-6.3.14). The controls were computed in a closed loop fashion, and 
the control computation to real-time ratio was found to be approximately 1:30. 

6 . 2 . 1 . 2 h-Boundary Layer Computations 

In order to obtain the pseudocontrol y* and the control a* in the h-boundary 
layer, we have implemented the successive approximati on scheme suggested in 
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Subsection 5.6.1. Specifically, we have performed the following iteration: 

(2.1) ap = 0 '' 


(2.2) Yi = 


COS 


-1 


1 + (T- Dfi-nL a? i) 

E mg 0 a 


(2.3) a^. = 


> i = UlO 


L +T 

a 


(2.4) a^. = - 


L aj + T sina'. - mg cosy. 
a 1 1 ^ ' 1 


L + T cosal 

CL 1 


We take Y* = Y2 q snd a* = a^g. Equation (2.2) is a complete form of (5.6.8), while 
(2.3) and (2.4) consist of a Newton iteration to solve for a in (5.6.2). If 
h>h*(E) we negate y* in order to achieve a dive. 

Figure 6.2.13 shows a plot if h*(E) versus E. This was obtained from the 
E boundary layer simulation of this subsection. We have also plotted, in Fig- 
ure 6.2.14, Y* 3S a function of h-h*(E) at two energy levels: E= 5482.7 m and 

E= 13258.9 m. Figures 6.2.13 and 6.2.14 represent the feedback control laws for 
the E and y boundary layers, respectively. 

6. 2. 1.3 Addition of a Predictive Feedback Term 

It was found that the system behaved better numerically when the nonlinear 
feedback terms described above were supplemented by a predictive term K^(h*-h). 
The control law is tany* = tany^^ + K^(h* - h) and a* satisfies (5.6.2). Addition 
of this predictive term has the same effect as imposing a penalty cost in addi- 
tion to the other cost criteria to penalize h when it is far from the "optimal" 
value of h computed as a control in the E-layer. Best system performance was 
obtained with = 0.005. 

6.2. 1.4 y-Boundary Layer Computations 

To obtain oc* for the y- boundary layer which is the lowest level of the hier- 
archical feedback structure, we added a predictive feedback term of the type 
K^(y*-y) subject to saturation limits of 0° and 12°. The use of SPT control 
together with predictive feedback terms stabilizes the control scheme and 
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alleviates inaccuracy due to numerical error in the control computation. Best 
performance was obtained for K^,= 0.6/cosa (a=76° while turning). 

6. 2. 1.5 The Initial Turn 

In Section 5.5 where the energy-state solution was obtained, transitions 
in altitude and velocity were instantaneous at constant energy by the very nature 
of the approximation. Hence it was cost-effective to zoom to the max-turn locus, 
complete the turn on the locus and zoom back to the optimal energy climb path. 

With h and y dynamics taken into account it is no longer clear whether this is 
the optimum strategy. Hence, a simulation of the aircraft trajectory was run 
using the controls calculated for the y-layer and two turning strategies were 
implemented and compared: (a) turn on max-turn locus and (b) turn at initial 

states. Table 6.2.1 presents these results together with simulation results for 
the flight path recommended by flight manuals. Also included in this table are 
the results of control laws using different gains for the predictive feedback 
terms, as well as the energy state solution. The results in the table indicate 
that the total time to the cruise arc is fairly insensitive to variations in 
these gains • The energy-state approximation assumes instantaneous 

maneuverability in h and y states and therefore provides a lower bound to the 
optimal solution. The best control strategy was found to be (b) above (Case 9 
in Table 6.2.1) and was only 4% above this lower bound while current practice 
(Case 10) was 21% above the lower bound and 17% above the best SPT solution (b). 
The strategy (a) (Case 8) which involves a climb to the max-turn locus for turning 
was 9% greater than (b). Hence although turning is faster on the locus (2.2 sec 
less), the initial distance from the locus (about 3000 m altitude here) makes 
(a) less efficient than (b). 

Figures 6.2.15, 6.2.16, 6.2.17, 6.2.18 and 6.2.19 give a view of the first 
20 seconds of the trajectory ((b) above), 10 seconds of which were taken up in 
turning. The complete climb to cruise arc is shown in Figures 6.2.20-6.2.30 
with the E-layer and h-layer trajectories included wherever they provide a good 
comparison with the y-layer to indicate the purpose of adding boundary layers. 

In particular. Figure 6.2.21 shows the zoom dive in altitude (h*) demanded by 
the E-layer and the well-rounded descent in altitude achieved by the y-layer. 
Figure 6.2.30, which is an al ti tude-Mach profile of the trajectory, includes in 
addition to the E-reduced order solution, the current practice trajectory which 



is recommended by' flight manuals. 

6.2.2 Cruise Arc 

Upon reaching the end of the energy-climb path of Section 6 , 2.13 some minor 
adjustments in altitude, velocity and flight path angle are necessary to attain 
the cruise values specified in Subsection 5.2.1. As mentioned in Subsection 
6.2. 1.1, there is a jump in velocity from 585. 2 to 590.2 m/sec at In 

the example simulation, y is slightly positive (4.4°) upon reaching and 
therefore must be reduced to 0°. In order to achieve cruise values for all the 
states simultaneously, a backward integration was performed beginning at cruise 
values maintaining E=0 and a as large as possible so as to increase y: these 

controls are u=l and a for E=0. In the real-time case, this integration could 
be done off-line and stored as a velocity-flight path angle profile. Forward 
integration was begun with the nominal cruise values (V= 585.2 m/sec, y=4.4°, 
etc.) using a=0 to decrease y and u for E = 0 and the states monitored to de- 
termine the intersection point with the backward trajectory to yield the optimum 
switching point. At this point the backward simulation controls are implemented: 
u = 1 and a for E=0 to arrive at the exact cruise values. In the example simu- 
lation, it took 17.2 secs to achieve these exact values.’ The same range could 
have been covered in 17.0 secs at V= 590.2 m/sec and y=0°, hence the cost of 
these corrections is very small, i.e., 0.2 secs. 

6.2.3 Descent from Cruise Arc to the Terminal Aircraft States 

This portion of the trajectory is a short range (30- 40 Km) optimization 
problem. It has been shown in Section 4.4 that the SPT approximation breaks 
down as the terminal states are approached. We therefore look to the energy- 
state approximation of Section 5.5 to provide clues as to the optimum path to 
the interception point. This approximation neglects h and y dynamics and so 
allows free (instantaneous) zoom dives and climbs. Hence, we find that the 
energy-state reduced order solution requires descent from the cruise arc along 
the max-velocity constraint until the final energy is reached after which there 
is a zoom climb (constant energy) to the required terminal altitude. In contrast 
to this solution we also have what is recommended in operations manuals for 
current supersonic aircraft: descent from the cruise arc along the max-velocity 

constraint until the final altitude is reached after which level (constant 
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altitude) flight is suggested to the required final energy (or velocity). 

Simulation of the trajectory along the max-velocity constraint was achieved 
by obtaining the two controls: alpha (a) and throttle (u) so as to satisfy 

where V_^',(E) is as in Figure A5.1.2 and V' „{E) is the derivative of with 
respect to E. 

V/e now consider a simulation of the zoom climb at constant energy from the 
max-velocity constraint to the final altitude as required by the energy-state 
reduced order solution. A period of pull-out from descent on the max-velocity 
constraint to ascent on the constant energy (E^) contour is necessary. Basically 
this involves changing the flight path angle from negative to positive as quickly 
as possible. Since at the final altitude, y^ = 0°, a switching strategy exists 
so that on the E^ contour, a maximum y is just reached after which the controls 
are switched so as to decrease y as fast as possible (d=0°) to satisfy the ter- 
minal constraint on y. There is also an optimum energy at which to begin pulling 
out of the descent so as to arrive onto the E^ contour with a nonnegative 
flight path angle. 

With a=0 on the energy contour (E^), it can be easily shown that 

(2.6) h = -E^ tan^y + h^ sec^y 

Hence, the switching height h^^ is given by 

(2.7) hs„ = -Ef 

Hence h^yj is the altitude corresponding to therefore a comparison 

of current altitude (h) to hc,i indicates whether has been reached. When 

h = hsw> ci must be set to zero so as to just bring y to zero at h^. The real- 
time algorithm is outlined below: 

(i) descend on constraint until pull-out energy with a and u from 
(2.5) 

(ii) pull-out with ct=cig, u=l until E = E^ 

(iii) maintain E=0 and allow y to increase until determined by 
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testing current altitude against 

(iv) set a=0 and u to maintain E=0 so that is attained at constant 
, energy (E^) with y=0 at 

(v) complete final turn with h=0, y=0, E=0 maximizing B. 

Example simulations were conducted using the above, algorithm to descend ; - 
from cruise to the final states. Two examples were implemented: 

(1) Ef = 1.2802X 10^ m, h^= 9. 144 x lo^ m, Yf=0°, 6f = 90°, M^ = 0.89 

(2) E^= 1.0668 X lo'^ m, = 9. 144 x lo^ m. Yf=0°, 6^ = 90°, = 0.56 

Note that horizontal plane constraints (i.e. , x.jr, y^) are satisfied by choosing 
the corredt cruise heading ($ ) and range on cruise (R ). The results for (1) 
are presented in Table 6.2.2 which also contains the energy-state reduced order 
solution and an Implementation of current practice. The effect of earlier pull- 
out from descent on the constraint is shown by the superiority of Case 4 

maX 

over Case 3 in Table 6.2.2. Note that as before the energy-state approximation 
solution provides a lower bound to the optimal solution since instantaneous 
maneuverability is assumed in the h and y states. The best real-time solution 
(Case 4) was 24% above this lower bound while current practice was 28% above 
the lower bound and about 4% above Case 4. 

The closeness of the result of Case 2 (current practice) to the best real- 
time solution (Case 4) motivated the simulation of example (2) above, the results 
of which appear in Table 6.2.3. Again the energy-state reduced order and current 
practice solutions are included. Here the best real-time solution (Case 4) is 
about 19% above the lower bound while current practice (Case 2) is 65% above 
the lower bound and about 38% above the best real-time solution. Hence taking 
the results of (1) and (2) together, current practice is almost as good as our 
best real-time solution in (1) where the final Mach number is high (M^ = 0.89) 
while in example (2) current practice is considerably less satisfactory than the 
best real-time solution. In (2) the final Mach number required is lower (M^ = 0.56) 
and we see that level flight through^ lower Mach number regions is time-consuming. 
Hence the optimum strategy for descent from the cruise arc is dependent upon 
the ppirticular final states required of the airplane. 

The whole trajectory of example (2) above, i.e.. Case 4 of Table 6.2.3, 
is displayed in Figures 6.2.31-6.2.43, of which Figures 6.2. 39 - 6.2.41 represent 
the three control (a, u, a) time-histories. We have added n'n the E reduced-order 
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solution and current practice wherever the comparisons are meaningful. Note 
that for the purpose of comparing the different cases in Table 6.2.3 a common 
range (68.86 Km) had to be considered so Case 4 includes extra range on cruise. 
With the exception of Figures 6.2.32, 6.2.42 and 6.2.43, this extra range is not 
shown in the figures. Common range considerations were also made in Tables 6.2.1 
and 6.2.2 with extra range added onto cruise wherever necessary. 

6*3 Numerica l Results Using a Real-Time Algorithm Based on the Energy State 
Approximation ^ -- - 

A real-time trajectory was obtained using the algorithm described in Sub- 
section 5.5.2 with turning only on the max-turn locus. Initial and final con- 
ditions chosen were: 


hQ = 3.353 X 10^ m 

h^ = 8.230X 10 

Mq = 0.6 


Eq = 5.517 X 10^ m 

E^ = 1.158 X 10' 

^0 " 

Xf = 212.0 Km 

o 

II 

o 

y^ = 169.6 Km 

Bq - 30° 

= 120° 


Figure 6.3.1 gives the energy changes with time along the trajectory. 
Roughly two thirds of the trajectory consists of reaching the cruise arc while 
the section from the cruise to the final conditions is by comparison much 
shorter: the irate of energy change being higher along the bank angle chatter 

arc. 

Figures 6.3.2 to 6.3.4 reveal the zoom dives and climbs which reflect the 
instantaneous changes in altitude and velocity allowed by the energy-state 
approximation. The real-time approximation used here which allows turniing only 
on the max-turn locus is clearly visible in Figure 6.3.5 which shows how the 
required 90° turn is broken up into an 8.3° initial turn and 81.7° final turn., 
Figure 6.3.6 shows the many variations in bank angle along the trajectory.^ By 
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the approximation mentioned above turning is with = and a = 0 for straight 


flight while 0 = ±min(a *a ) during chattering to maintain average 6 = 0. 

Ill d 


The 


thrust control history is quite simple: 


T = all the way to the cruise arc. 


T,= T. < while cruising since the cruise arc is an intermediate thrust singular 

C maX 

arc followed by T=0 while the bank angle chatters and then T = i^or the final 
turn. , 


The horizontal plane projections of the trajectory appear in Figures 6.3;8" 
6.3.11. Figure 6.3.11 indicates that a very small portion of the trajectory, 
time is spent in turning on the max-turn locus as compared to the straight- flight 
portions. The alpha control history is given in Figure 6.3.12, The altitude- 
Mach number profile of Figure 6.3.13 provides a very interesting view of the 
trajectory. From the initial conditions which are on the max-turn locus, the 
trajectory executes the initial turn followed by a zoom dive to the minimum time 
energy climb path which it follows to the cruise point. After cruising the bank 
angle is allowed to chatter while moving down the max velocity constraint. 

The trajectory leaves the constraint by a zoom climb to the max-turn locus where 
it executes the final turn to reach the end conditions which are also on the 
locus. If the initial or final conditions were not on the locus, the trajectory 
would include zoom dives or climbs to reach the locus. Figure 6.3.14 is an 
altitude-range profile where the rapid changes in altitude are clearly visible 
while Figure 6.2.15 shows the changes in Mach number with range. 

In order to perform the t^ iterations of the algorithm presented in Sub- 
section 5.5.2, some assumptions have to be made about the target dynamics. We 
assume that it moves in a straight line at constant altitude and velocity. Then 
the steps in the algorithm are as follows: 

(1) obtain initial estimate of interception point by calculating the 
shortest distance to target path from pursuer 

(2) evaluate t^=time for target to reach this point 

(3) implement (2) - (12) of algorithm in Subsection 5.5.2 

(4) compute pursuer trajectory time (t^) to interception point 

(5) define the difference as = and perform approximate Newton 

iteration as follows: 



f(t 


f ^ 
^old 


old 


y " ^1 
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(6) go to (3) if accuracy on is not met. 

(1) - (6) above were implemented for two cases given in Figures 6.3.16 and 
6.3.17. The first case required six iterations for a 180° turn to an accuracy 
of 0.1 secs. The second case (Figure 6.3.17) required 11 iterations. This was 
the worst’ case: 90° turn to a target which has already passed the point where 

its path is at the shortest distance from pursuer. Hence the two cases indicate 
good convergence. Note that a better initial estimate of the interception point 
(step (1) above) using proportional navigation or some other scheme would reduce 
the number of iterations. 
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Table 6.2.1 Summary of Simulation Results for Trajectories to Cruise Arc 


CASE 

Time to cruise arc, 
starting at H = 593.4 m, 

M = 0.9, energy = 5482.7 m 
Final range =182.9 Km 
altitude = 12192 m, 

Mach = 2.0, energy = 29949.3 m 

REMARKS 

1. Energy-layer 

380.9 secs 

This approximation treats 
h and y as controls and 
holds the derivatives Fi 
and Y equal to 0, In- 
tuitively, it assumes that 
h and y change much more 
quickly than E. That is, 
h and y are characterized 
by sudden changes (boundary 
layers) for short times 
and almost steady behavior 
for the rest of the time. 

2. h-layer with 
nonlinear and 
linear feedback 
gain on tany - 
1 inear gain = 0.01 
on (H*(E)-H). 

387.1 secs 

This approximation treats 
Y as a control, assuming 
that Y varies faster than 
h does. Linear feedback 
is also used. It has the 
same effect as imposing a 
quadratic penalty cost in 
addition to the other cost 
criteria to penalize h when 
it is far from the optimal 
h-control of the E-layer. 

3. h-layer with 

386.9 secs 

See remarks above. 

4. Y- layer with 
^GAIN ~ 

linear feedback 
gain on a, 

^GAIN"°'^ 

0<a < 12° 

397.2 secs 

The control uses a nonlinear 
SPT control together with 
linear feedback terms which 
represent quadratic penalties 
for being away from "optimal" 
values of h computed as a 
control in the E- layer and 
of Y computed as a control 
in the h layer. This pro- 
cedure stabilizes the control 
scheme and alleviates in- 
accuracy due to numerical 
error in the control com- 
putation. 
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Table b.Z . 1 (cont.) 


CASE 

TIME TO CRUISE ARC 

. REMARKS 

5. Y-laye>" with 
0 <a< 12” 

397.1 secs 

See remarks in 4. 

6. Y“ layer 
0<a< 12” 

398.1 secs 

See remarks in 4. 

7, Y“ layer 
^GAIN 

^GAIN^®*^ 

0 <a < 12” 

397.0 secs 

See remarks in 4. 

8. Y' layer with 90° 
turn on max-turn 
locus 

0<a < 12° 

a = 76° (on locus) 

438.5 secs 

Path to locus simulated 
using a flip-flop control 
and Y constrained to be less 
than 47°. Aircraft flies to 
the maximum rate turn locus 
to turn. 

9. Y-layer with 90° 
turn begun at 
initial state 
Yq = 0.6/cosa 

hg = 0.005 

0<a< 12° 

a = 76° while turninc 

403.0 secs 

The aircraft executes its 
turn without returning to 
the max rate turn locus. 

This strategy is much more 
efficient than (8) given the 
initial distance from the 
locus (about 3000 m altitude). 

. . 

10. Current practice: 
turn using cr~ 60° , 
at constant velo- 
city, constant 
altitude acc. to 
M=0.9, constant 
Mach § cl i mb to 
12192 m and 
Y = 21.8”, constant 
altitude acc. to 
cruise point. 

472.7 secs 

■ 

This scheme is typical of 
paths recommended in operations 
manuals for current super- 
sonic aircraft. 
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Table 6.2.1 (cont. ) 


CASE 

TIME TO CRUISE ARC 

■ 

REMARKS 

11. Energy layer with 

388.7 secs 

See remarks for 1. The 

zoom to max turn 


additional 7.8 secs is the 

locus, 90° turn 


time taken for a 90° turn 

on the locus and 
zoom back to 
energy layer 
{Parsons solution) 


on the max turn locus. 



Table 6.2.2 Sunmary of Simulation Results for Terminal 
Trajectories off the Cruise Arc 


■ . Case 

Time from cruise arc: 
energy=29946.6 m, h=12192 m 
M=2.0 to final states: 
energy- 12801. 6 m, M=0.89, 
h=9144 m, Yt0°,^. . , 

range=46.19 Kirt " ' 

Remarks 

1. Energy-state approxi^ 
mat! on solution: 
down max -velocity 
constraint to E^= 1280 1.6 
m; zoom cl imb to 
h=9144 m at constant 
energy; 90^ turn at 
E=12801.6 m, h=9144 m 

92.5 secs 

See Parsons (1972) 
for this solution. 
This is a lower 
bound for the optimal 
solution since in- 
stantaneous maneu- 
verability is assumed 
in h and y states 

2. Current practice: 
down max-veloci ty 
constraint until 
h=9144 m and y=0°; 
level flight at 
h=9144 m to E=12801 m, 
M=0.89; 90° turn with 
o, a, u sg as to 
maintain E=0, h=0, y=0 

118.6 secs 

This scheme is 
typical of paths 
recormiended in oper- 
ations manuals for 
current supersonic 
aircraft. 

3. Real-time solution; 
down max-velocity 
constraint till 
energy=12801.6 m; 
pull-out with a=a , 

(-21°)= 

Y=0 with u=l till 
energy= 1280 1.6 m again; 
E=0, a=0 so as to obtain 
Y=0 at h=9144 m; 90° 
turn with a, .a, u.so as 
Jo maintain E=0, h=0, 

Y=0 

130.9 secs 

This solution at- 
tempts to follow (1) 
as closely as pos- 
sible; Y and h dy- 
namics require a 
period of pull-out 
from descent on max- 
velocity constraint. 
The final state Yf=0° 
requires y to be ‘ 
limited to 
while on 
energy contour. 

4. Best real-time solution: 
down max-velocity 
constraint till 
energy=15849.6 m; pull- 
out with ct=o. t u=l till 
energy= 1280 176 m; main- 
tain E=0 till Y reaches 
■^max (“'ll")! 
with a=0 so as to reach 
h=9144 m with y=0°; 

90° turn with a, a, u 
§0 as to maintain E=0, 
h«0, Y=0 

114.6 secs 

This solution is 
better than (3) be- 
cause of earl ier 
pull-out from 
descent on V 
constraint 
(3048 m sooner). 
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Table 6.2.3 Summary of Simulation Results for Terminal 
Trajectories off the Cruise Arc ; 


Case 

Time from cruise arc: 
Energy=29946.6 m, h= 
12192 m, M=2.0 to final 
states: Energy= 10668 m 

M=0.56, h=9144 m, 

Y=0°, range=68.86 Km 

Remarks 

1. Energy-state approximation 
solution: down max-velocity 

constraint to E^= 10668 m; 
zoom climb to h=9144 m 
at constant energy; 90° 
turn at E=10668 m, h= 

9144 m. 

138.8 secs 

See Parsons (1972) for 
this solution. This 
is a lower bound for 
the optimal solution 
since instantaneous 
maneuverability is 
assumed in h and y 
states. 

2. Current practice: dov^n max- 

velocity constraint until 
h=9144 m and y=0°; level 
flight at h=9144 m to E= 
10668 m, M=0.56; 90° turn 
with CT,a,u^so as to main- 
tain E=0, h=0, Y=0. 

228.0 secs 

This scheme is typical 
of paths recommended 
in operations manuals 
for current supersonic 
aircraft. 

3. Real-time solution: down 

max-velocity constraint 
till E=13716 m; pull-out 
with u=l till^E= 

10668 m, maintain E=0 
till Y. reaches Y^a. (=59.6°); 
keep E=0 with a-D so as to 
reach h=9144 m with y=0°; 

90° turn with g,ct,u^so as to 
maintain £=0, h=0, y=0 

169.1 secs 

This solution follows 
(1) with a pull-out 
from descent on V 
constraint 10000 
ft before reaching 
E^=10668 m. 

4. Best real-time solution; 
down max-velocity constraint 
till E=15240 m; (earlier 
pull-out than in (3)); pull- 
out with a=a , u=l till E= 
10668 m; maintain £=0 till 
Y reaches y„,,„ (=59.6°); 
keep t=0 so ^as to reach 
h=9144 m with y=0°; 90° 
turn with a,a,u so as to 
maintain £=0, h=0, y=0. 

165.2 secs 

This solution is better 
than (3) because of 
earlier pull-out (1524 
m sooner than (3)) 
from descent on 
constraint 
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Figure 6.2.21 Altitude Time-History. h*(E) denotes the desired altitude 
history from the E-boundary layer. 
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Figure 6.2*26 v-Time Histories. Y*(h) denotes the desired ^-history computed 
from the h-boundary layer. 
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Figure 6.2.27 a-Control Histories. a,{E) is the control history computed by the 

energy boundary layer. * o! 2 (h) is the control history computed by the 
h-boundary layer. The actual control used by the aircraft is a. 
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Figure 6,3.7 Changes in Thrust with Time 
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Figure 6.3.10 Range in Horizontal Plane with Time 
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Figure 6.3.16 Convergence of t^ Iterations (180° turn) 
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Figure 6.3.17 Convergence of t^: Iterations (90° turn) 



6.4 S umma ry o f ]tae r i c a 1_ Si mul ation Results 


Table 6.4.1 Case 1 - Full Trajectory Simulation 
Using SPT Algorithms 


In itial States 
Eq = 5.483 X 103 m 
hQ = 0.594 X 1Q3 m 
- 0.92 

To = 0° 

Xq = 0 Km 
yg = 0 Km 
So = 0» 


Fin al States 

= 1.067 X lo'^ m 
= 9.144 X 1Q3 m 
= 0.57 

Tf = 0° 

= -2.03 Km 
y^ = 277.8 Km 
= ISO'^ 


h 

/h 


5 i mula tion Results 

The six state simulation was done in 3 parts as below with total 
(minimum) time =614.0 secs using the SPT algorithms of Section 6.2. 


Part 

Optimal Control Law 

Minimum Time (sec) 

I. Climb to 
cruise arc 

Case 9 of Table 6.2.1. 
Ax = 0.43 Km, 

Ay = 180.7 Km 
Figures 6.2.15-6.2.30 

399.3 

II. Minor 

Corrections 
to achieve 
cruise values 
and cruise 
portion 

See Section 6.2.2. 
Ax = 0 , Ay = 64 . 4 Km 

109.3 

III . Descent 

from cruise 

Case 4 of Table 6.2.3 
Ax =-2.46 Km, 


arc 

Ay = 32.7 Km 

Figures 6.2.31 - 6.2.43 

105.4 



, T^b.le;,,6.4.2 v:£^se.J2. - :,G41mb to Cruise Are Simulation.. 
'f'^-Us'ing 5 P[f Algorithm . V 


Initial States ; . 
Same as in Case 1 
(T^bie 6.4a),, 


Final /States- : ; 

- 2.995 x- io'+ .m 
hf = 1.219.x 10"^ .m 

= 2.00 

Jf = 0 ° 

- 0.43 Km . . 
y.p = 190.78 Km . 
3.p = 90° 


Simulation Results 


This is a six state simulation and consists of Part 1 of Table 6.4.1 
plus "the minor corrections to achieve cruise values" of Part .2 of Table 
6.4.1 which took 17.2 secs and covered 10.04 Km. Hence total (minimum) 
time = 416.5 secs. The SPT algorithm of Section 6.2 was used. 


3.1|6 



Table 6.4.3 Case 
Using 

3 - Descent from Cruise Arc Simulation 
SPT Algorithm (M.p = 0.56) 

Initial States 

Final States 

Eg = 2.995 X lo** m 

= 1.067 X lo** m 

hg = 1.219 X 10“ m 

h.p * 9.144 X 10^ m 

Mg = 2.00 

= 0.56 

O 

O 

If 

O 

yf = 0® 

Xq = 0 Km 

x^ = 32.7 Km 

yg = 0 Km 

= 2.46 Km 

3g - 0» 

= 90*=^ 


S imulation Results 

This six state simulation consists of Part 3 of Table 6.4.1 taking 
105.4 secs and uses the SPT algorithm of Section 6.2. 


Table 6.4*4 Case 4 - Descent from Cruise Arc Simulation 
Using SPT Algorithm =0.89) 


Initial States 
Same as in 
Case 3 

(Table 6.4.3) 


Final States 

= 1.280 X lO'^ m 

hn = 9.144 X 10^ m 
t 

= 0.89 

Yf = 0° 

= 37.4 Km 
= 2.46 Km 
= 90° 


Simulation Results 

This six state simulation uses the same algorithm used in Table 6.4.3 
and is described in detail as Case 4 of Table 6.2.2. Total (minimum) time 
is 99.85 secs for a range of 37.5 Km. 
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Table 6.4.5 Case 5 - Target Interception (180° turn) 
Using the Energy-State Approximation 


Pursuer States 


Eq = 6.706 X 1q 3 m 
Hq = 3.176X 103 m 

Mq = 0.8 

YO = 0» 

Xq = 0 Km 
Yq = 0 Km 


Target States 


Eq = 1.600 X 10^+ m 
Mq = 9.144 X 103 m 
Mg = 1.21 

Y„ = 0° 


Xq = 212.0 Km 
= 212.0 Km 


3q = 180‘ 


= 1.707 X 10^ m 
h^ = 9 . 144 X io3 m 
= 1.3 

Yf = 0° 

x.^ = 46.6 Km 
y^ = 212.0 Km 
= 180*^ 


E^ = 1.600 X 10^ m 
= 9 . 144 X lo3 m 
= 1.21 
= 0 ° 


X.J, = 46.6 Km 
y^ = 212.0 Km 

3f = 180° 


Simulation Results 

This simulation uses the energy-state approximation 
(4 state model of Section 5.5). The algorithm of Sec- 
tion 6.3 is used to iterate on tf which is the time to 
interception. Six iterations were required (Figure 
6.3.16) for convergence to an accuracy of 0.1 secs to 
yield a minimum-time trajectory taking 452.1 secs to 
interception. The trajectory reflects the energy-state 
approximation by the zoom maneuvers in altitude and velocity--the whole trajec- 
tory IS similar in character to that represented by Figures 6.3.1 -6 3 15 
which used different initial and final conditions from those stated above 
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Table 6.4.6 Case 6 - Target Interception (90° turn) 
Using the Energy-State Approximation 


Pursuer States 


Eq = 6.706 X 1q 3 m 
hQ = 3.176 X iq 3 m 

Mq = 0.8 

Yo = 0° 

Xq = 0 Km 
Vn = 0 Km 


Ef = 1.707 X 10^* m 
hf = 9.144 X iq 3 m 


= 1.3 


Xf = 212.0 Km 
yf = 294.6 Km 
3f = 90° 


Target States 


Eq = 1.600 X lO"* m 
hQ = 9.144 X io3 m 
Mq = 1.21 

Y(j = 0» 

Xq = 212.0 Km 
= 42.4 Km 


Ef = 1.600 X 10"+ m 
= 9.144 X iq 3 m 
= 1.21 
Yf = 0° 

Xf = 212.0 Km 
y^ = 294.8 Km 


3q = 90° 


Simulation Results 


3f = 90< 


Taj^g^t 


This simulation is similar to that of Table 
6.4.5 in that the energy-state approximation 

is used and the same algorithm (Section 6.3) p >> 

is also used but differs in the initial and ^ 

final states. 11 iterations were required (Figure 6.3.17) and the minimum- 
time trajectory takes 690.3 secs to interception. 
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CHAPTER 7 
CONCLUSIONS 


This report has described work performed on the development of a hierarchi- 
cal real-time algorithm for optimal three-dimensional aircraft maneuvers using 
Singular Perturbation Theory (SPT). New theoretical results justify and develop 
systematic methods for real-time computation of nonlinear feedback controls, by 
means of SPT and provide an assessment of the accuracy of the resulting SPT con- 
trol. Practical results apply SPT to obtain a real-time feedback law for the 
three-dimensional minimum time long range intercept problem for an F-4 aircraft 
model (six state, three control variable, point mass model). Nonlinear feed- 
back lav/s are presented for computing the optimal control variables u (throttle), 
a (bank angle) and a (angle-of-attack) as a function of target and pursuer air- 
craft states and desired terminal conditions. The SPT control law results in a 
hierarchical nonlinear feedback structure. It is supplemented by predictive 
feedback terms for small deviations from the optimal trajectory and for maneuvers 
near the terminal time where the SPT approximation is not valid. The F-4 simu- 
lation results using the SPT control law show minor sacrifice in accuracy over 
the off-line optimization results, in the long range intercept case. 

A real-time capability assessment of the SPT algorithm on a microcomputer 
has been performed and based on the results presented in this report, it may be 
concluded that real-time, three-dimensional long range aircraft trajectory 
optimization is possible using SPT, The implementation of this algorithm on a 
microcomputer is estimated to result in a control update cycle time of 27 msec, 
which is almost four times smaller than the common radar sampling interval of 
100 msec. The storage and computational requirements of the algorithm are found 
to be well suited for on-board real-time implementation on a microcomputer. 

The accuracy of the SPT solution is analyzed and it is shown how "continua- 
tion-type" methods may be used to obtain exact optimal trajectories starting 
from the SPT solution. The advantage of using predictive terms to supplement 
the SPT feedback laws is demonstrated for the aircraft trajectory optimization 
problem. In particular, it is shown that the SPT approximation breaks down near 
the terminal phase and must be corrected by "continuation" and Generalized 
Multiple Scale (GMS) methods. 
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